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1 MATH FACTS 

1.1 Vectors 

1.1.1 Definition 

We use the overhead arrow to denote a column vector, i.e., a number with a 
direction. For example, in three-space, we write 




The elements of a vector have a graphical interpretation, which is particularly 
easy to see in two or three dimensions. 



z 




1. Vector addition is pointwise. 



a + b = c 



f 2 1 




f 2 1 




f 5 1 


1 


► + < 


3 


> = < 


4 


1 7 J 




. 2 . 




i 9 j 



Graphically, addition is stringing the vectors together head to tail. 
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2. Scalar multiplication is pointwise. 




1.1.2 Vector Magnitude 

The total length of a vector of dimension m, its Euclidean norm, is given by 



m 




this scalar is commonly used to normalize a vector to length one. 



1.1.3 Vector Dot Product 

The dot product of two vectors is the sum of the products of the elements: 



x-y = x^y = 

i=l 

The dot product also satishes 

X ■ y = \ \x\\\\y\\ cos 9, 

where 6 is the angle between the vectors. 

1.1.4 Vector Cross Product 

The cross product of two three-dimensional vectors is another vector, xxy = z, 
whose 

1. direction is normal to the plane formed by the two vectors, 

2. direction is given by the right-hand rule, rotating from x to y, 




1.2 Matrices 
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3. magnitude is the area of the parallelogram formed by the two vectors - 
the cross product of two parallel vectors is zero - and 

4. (signed) magnitude is equal to | |a;| 1 1 \y\ \ sin 9, where 9 is the angle between 
the two vectors, measured from x to y. 

The schoolbook formula is 

f X2y^ - x^y2 1 
xxy = I xsyi - xiya > . 

[ xm - X2yi ] 



1.2 Matrices 

1.2.1 Definition 

A matrix, or array, is equivalent to a set of row vectors, arranged side by side, 
say 



A = 




2 3 
1 3 
7 2 



This matrix has three rows (m = 3) and two columns (n = 2); a vector is 
a special case of a matrix with one column. Matrices, like vectors, permit 
pointwise addition and scalar multiplication. We usually use an upper-case 
symbol to denote a matrix. 



1.2.2 Multiplying a Vector by a Matrix 

If Aij denotes the element of matrix A in the i’th row and the j’th column, 
then the multiplication c = Av is constructed as: 



Cj — AiiVi + Ai2V2 + ■ ■ ■ + AinVn — ^ 

i=i 



where n is the number of columns in A. c will have as many columns as 
A has rows (m). Note that this multiplication is well-dehned only if v has 
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as many rows as A has columns; they have consistent inner dimension n. 
The product vA would be well-posed only if A had one row, and the proper 
number of columns. There is another important interpretation of this vector 
multiplication: Let the subscript ; indicate all rows, so that each A.j is the 
j’th column vector. Then 



C = Av = A:iVi A.,2V2 H h A:nVn- 

We are multiplying column vectors of A by the scalar elements of v. 

1.2.3 Multiplying a Matrix by a Matrix 

The multiplication C = AB is equivalent to a side-by-side arrangement of 
column vectors = AB.j, so that 



C = AB= [AB,i AB,2 ■ ■ ■ AB,k], 

where k is the number of columns in matrix B. The same inner dimension 
condition applies as noted above: the number of columns in A must equal the 
number of rows in B. Matrix multiplication is: 

1. Associative. {AB)C = A{BC). 

2. Distributive. A{B + C) = AB + AC, {B + C)A = BA + CA. 

3. NOT Commutative. AB ^ BA, except in special cases. 

1.2.4 Common Matrices 

Identity . The identity matrix is usually denoted I, and comprises a square 
matrix with ones on the diagonal, and zeros elsewhere, e.g.. 



-^3x3 — 



1 0 0 
0 1 0 
0 0 1 



The identity always satishes AJ„xn = Imxm^ = A. 




1.2 Matrices 
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Diagonal Matrices . A diagonal matrix is square, and has all zeros off the 
diagonal. For instance, the following is a diagonal matrix: 



A 



4 0 0 

0-2 0 
0 0 3 



The product of a diagonal matrix with another diagonal matrix is diagonal, 
and in this case the operation is commutative. 



1.2.5 Transpose 

The transpose of a vector or matrix, indicated by a T superscript results 
from simply swapping the row-column indices of each entry; it is equivalent to 
“flipping” the vector or matrix around the diagonal line. For example. 



I 2 I — = {12 3} 



■ 1 2 ■ 






4 5 


A^ = 


14 8 
2 5 9 


8 9 







A very useful property of the transpose is 

(ABf = B^A^. 



1.2.6 Determinant 

The determinant of a square matrix A is a scalar equal to the volume of the 
parallelepiped enclosed by the constituent vectors. The two-dimensional case 
is particularly easy to remember, and illustrates the principle of volume: 



det(A) — A 11 A 22 ~ A 21 A 12 




det 



1 + 1 = 2 . 
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In higher dimensions, the determinant is more complicated to compute. The 
general formula allows one to pick a row k, perhaps the one containing the 
most zeros, and apply 



j=n 

det{A) = J2Akj{-l)’^+^Ak„ 

i=i 

where A^j is the determinant of the sub-matrix formed by neglecting the fc’th 
row and the j’th column. The formula is symmetric, in the sense that one 
could also target the fc’th column: 

j=n 

det{A) = Y^A^,{-l)’^+^A,k. 

i=i 

If the determinant of a matrix is zero, then the matrix is said to be singular - 
there is no volume, and this results from the fact that the constituent vectors 
do not span the matrix dimension. For instance, in two dimensions, a singular 
matrix has the vectors colinear; in three dimensions, a singular matrix has 
all its vectors lying in a (two-dimensional) plane. Note also that det{A) = 
det{A^). If det{A) ^ 0, then the matrix is said to be nonsingular. 

1.2.7 Inverse 

The inverse of a square matrix A, denoted A~^, satisfies AA~^ = A~^A = I. 
Its computation requires the determinant above, and the following dehnition 
of the n X n adjoint matrix: 



ad] {A) 



(-l)'+'An 






(-l)'+"Ai„ 



n T 



{-iT+^Ann 



1.2 Matrices 
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Once this computation is made, the inverse follows from 

^-1 ^ adjjA) 
det{A) 

If A is singular, i.e., det{A) = 0, then the inverse does not exist. The inverse 
hnds common application in solving systems of linear equations such as 

Ax = b — X = A~^b. 



1.2.8 Trace 

The trace of a matrix is simply the sum of the diagonals: 



tr{A) = '^Aii. 

i=l 

1.2.9 Eigenvalues and Eigenvectors 

A typical eigenvalue problem is stated as 



Ax = AT, 

where A is an n x n matrix, T is a column vector with n elements, and A is 
a scalar. We ask for what nonzero vectors x (right eigenvectors), and scalars 
A (eigenvalues) will the equation be satished. Since the above is equivalent to 
(A — \I)x = 0, it is clear that det{A — XI) = 0. This observation leads to the 
solutions for A; here is an example for the two-dimensional case: 



A 

A- AJ 
det{A — XI) 



4-5 ^ 

2 -3 J ' 

■ 4 - A -5 ■ _ 

2 -3- A 

(4-A)(-3-A) + 10 
A^ - A - 2 
(A -h 1)(A - 2). 
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Thus, A has two eigenvalues, Ai = —1 and A 2 = 2. Each is associated with a 
right eigenvector x. In this example. 



' 5 
2 



- Ai/)fi 




Xi 



0 — ^ 

0 — ^ 

{^ 2 / 2 , V 2 / 2 Y . 



{A 
' 2 
2 



— \2I)X2 




X2 



0 — ^ 

0 — ^ 

, 2v^/29}^. 




Eigenvectors have arbitrary magnitude and sign; they are often normalized to 
have nnity magnitnde, and positive first element (as above). A set of n eigen- 
vectors is always linearly independent. The condition that rank{A — Xil) = 
rank{A) — 1 indicates that there is only one eigenvector for the eigenvalue A*. If 
the left-hand side is less than this, then there are mnltiple nniqne eigenvectors 
that go with A*. 

The above discussion relates only the right eigenvectors, generated from the 
eqnation Ax = Xx. Left eigenvectors, also useful for many problems, pertain 
to the transpose of A: A^y = Xy. A and share the same eigenvalues A, 
since they share the same determinant. Example: 



{A^ 

5 

-5 



- Ai)yi 




Vi 



0 — ^ 

0 — ^ 

-5v^/29}^. 



{2v^/29, 



(A^ 

2 

-5 -5 



X2l)y2 
2 

■ V2 



V2 



0 ^ 

0 — ^ 

{^ 2 / 2 , -V 2 / 2 Y . 
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1.2.10 Modal Decomposition 

The right and left eigenvectors of a particular eigenvalue have unity dot prod- 
uct, that is xfiji = 1, with the normalization noted above. The dot product 
of a left eigenvector with the right eigenvector of a different eigenvalue is zero. 
Thus, if 

X = [xi - ■ ■ T„] , and 
Y = 

then we have 

Y^X = /, or 
= X-\ 

Next, construct a diagonal matrix of eigenvalues: 

■ Ai 0 

A = 

0 \n 

it follows that 

AV = 1/A - 
A = VAW^ 

n 

1=1 

Hence A can be written as a sum of modal components.^ 

1.2.11 Singular Value 

Let G{s) be an m X n, possibly complex matrix. The singular value decompo- 
sition (SVD) computes three matrices satisfying 

^By carrying out successive multiplications, it can be shown that has its eigenvalues 
at and keeps the same eigenvectors as A. 
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G = UTV*, 

where U is m x m, T is m x n, and V is nx n. The star notation indicates a 
complex-conjugate transpose. The matrix S is diagonal, with the form 

’ (Ti 0 0 0 ’ 

0-00 
^ “ 0 0 dp 0 ’ 

.0000. 

where p = min{m,n). Each nonzero entry on the diagonal is a real, positive 
singular value, ordered such that ai > a2 > ■ ■ ■ Cp. The notation is common 
that (Ti = a, the maximum singular value, and dp = a, the minimum singular 
value. The auxiliary matrices U and V are unitary, i.e., they satisfy X* = 
X~^. Like eigenvalues, the singular values of G are related to projections, cij 
represents the Euclidean size of the matrix G along the dth singular vector; 



a = maa;||a;||=i| |Ga; 

a = mm||a,||=i||Ga;| 

Other properties of the singular value include: 

• a{AB) < a{A)a{B). 

• t{A) = ^J\^ax{,A*A). 

• Q_{A) = ^J\jnin{,A*A). 

• Qi{A) = l/a{A~^). 

• T{A) = l/g_{A~^). 




1.3 Laplace Transform 
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1.3 Laplace Transform 

1.3.1 Definition 

The Laplace transform converts time-domain signals into a frequency-domain 
equivalent. The signal y{t) has transform T(s) defined as follows: 



fOO 

F(s) = L{y{t)) = / y{r)e-^^dT, 

Jo 

where s is an unspecified complex number; Y (s) is considered to be complex as 
a result. Note that the Laplace transform is linear, and so it is is distributive: 
L{x(t) + y(t)) = L{x(t)) + L{y(t)) will hold throughout. The following table 
gives a list of some useful transform pairs and other properties, for reference. 

The last two properties are of special importance: for control system design, 
the differentiation of a signal is equivalent to multiplication of its Laplace 
transform by s; integration of a signal is equivalent to division by s. The other 
terms that arise will cancel if ?/(0) = 0, or if ?/(0) is finite, so they are usually 
ignored. 

1.3.2 Convergence 

We note first that the value of s affects the convergence of the integral. For 
instance, if y{t) = e*, then the integral converges only for Re{s) > 1, since 
the integrand is in this case. Convergence issues are not a problem for 
evaluation of the Laplace transform, however, because of analytic continuation. 
This result from complex analysis holds that if two complex functions are 
equal on some arc (or line) in the complex plane, then they are equivalent 
everywhere. This fact allows us to always pick a value of s for which the 
integral above converges, and then by extension infer the existence of the 
general transform. 

1.3.3 Convolution Theorem 

One of the main points of the Laplace transform is the ease of dealing with 
dynamic systems. As with the Fourier transform, the convolution of two sig- 
nals in the time domain corresponds with the multiplication of signals in the 
frequency domain. Consider a system whose impulse response is g{t), being 
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1.3 Laplace Transform 
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driven by an input signal x(t); the output is y(t) = g(t)*x(t). The Convolution 
Theorem is 



y{t) = f g{t -T)x{r)dT = G{s)X{s). 

JO 

Here’s the proof given by Siebert: 



roo 

F(s) = / y(t)e~^^dt 

Jo 

roo r rt 

= / 9{t-T) x{t) dr 

Jo Uo J 

roo r roo 

= / g{t - r) l{t - t) x{r) dr 

Jo Uo 

roo 



dt 



e ^^dt 



roo roo 

/ x{t) / g(t — t) — t) dt 
Jo L^o 



dr 



roo 

= / x{r) G{s)e~^'^ dr 

= G{s)X{s) 

When g{t) is the impulse response of a dynamic system, then y{t) represents 
the output of this system when it is driven by the external signal x{t). 



1.3.4 Solution of Differential Equations by Laplace Transform 

The Convolution Theorem allows one to solve (linear time-invariant) differen- 
tial equations in the following way: 

1. Transform the system impulse response g{t) into G(s), and the input 
signal x{t) into W(s), using the transform pairs. 

2. Perform the multiplication in the Laplace domain to hnd T(s). 

3. Ignoring the effects of pure time delays, break Y (s) into partial fractions 
with no powers of s greater than 2 in the denominator. 

4. Generate the time-domain response from the simple transform pairs. 
Apply time delay as necessary. 

Specihc examples of this procedure are given in a later section on transfer 
functions. 
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2 KINEMATICS OF MOVING FRAMES 

2.1 Rotation of Reference Frames 

We say that a vector expressed in the inertial frame has coordinates x, and 
in a body-reference frame Xb- For the moment, we assume that the origins of 
these frames are coincident, but that the body frame has a different angular 
orientation. The angular orientation has several well-known descriptions, in- 
cluding the Euler angles and the Euler parameters (quaternions). The former 
method involves successive rotations about the principle axes, and has a solid 
link with the intuitive notions of roll, pitch, and yaw. Quaternions present a 
more elegant and robust method, but with more abstraction. We will develop 
the equations of motion using Euler angles. 

Tape three pencils together to form a right-handed three-dimensional coordi- 
nate system. Successively rotating the system about three of its own principle 
axes, it is easy to see that any possible orientation can be achieved. Eor ex- 
ample, consider the sequence of [yaw, pitch, roll]: starting from an orientation 
identical to some inertial frame, rotate the movable system about its yaw axis, 
then about the new pitch axis, then about the newer still roll axis. Needless 
to say, there are many valid Euler angle rotation sets possible to reach a given 
orientation; some of them might use the same axis twice. 




Eigure 1: Successive application of three Euler angles transforms the original 
coordinate frame into an arbitrary orientation. 



A Erst question is: what is the coordinate of a point hxed in inertial space, 
referenced to a rotated body frame? The transformation takes the form of 
a 3x3 matrix, which we now derive through successive rotations of the three 




2. 1 Rotation of Reference Frames 
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Euler angles. Before the first rotation, the body-referenced coordinate matches 
that of the inertial frame: = x. Now rotate the movable frame yaw axis 

{z) through an angle 0. We have 

cos 0 sin 0 0 

= — sin0 COS0 0 x^ = R{(j))x^. (1) 

[ 0 0 1 J 

Rotation about the 2 ;- axis does not change the 2 ;-coordinate of the point; the 
other axes are modihed according to basic trigonometry. Now apply the second 
rotation, pitch about the new y-axis by the angle 6: 

cos 9 0 — sin 0 

x^= 0 1 0 x,^ = R{9)x^\ (2) 

sin 6 0 cos 6 

Finally, rotate the body system an angle 0 about its newest x-axis: 

r 1 0 O' 

x^ = 0 COS0 sin 00 x^ = R{'ip)x^ . (3) 

0 — sin0 COS0 

This represents the location of the original point, in the fully-transformed 
body-reference frame, i.e., x^. We will use the notation x^ instead of x^ from 
here on. The three independent rotations can be cascaded through matrix 
multiplication (order matters!): 

Xb = R(0)i?(0)R(0)T (4) 

c9c(f) c9s(j) —s9 

= —c-ipstp + s-ipsOc^) c'lpccf) + s%Ijs6s( 1) stpcO x 
s'lps^) + cifjsOccI) —s-ipccj) + c'lpsOs'ip cipc9 

= R{4>,6,ip)x. 

All of the transformation matrices, including R(0, 6, 0), are orthonormal: their 
inverse is equivalent to their transpose. Additionally, we should note that the 
rotation matrix R is universal to all representations of orientation, including 
quaternions. The roles of the trigonometric functions, as written, are specihc 
to Euler angles, and to the order in which we performed the rotations. 

In the case that the movable (body) reference frame has a different origin than 
the inertial frame, we have 
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X = xq + R^Xb, (5) 

where xq is the location of the moving origin, expressed in inertial coordinates. 

2.2 Differential Rotations 

Now consider small rotations from one frame to another; using the small angle 
assumption to ignore higher-order terms gives 



1 


S(j) 


-S9 


—S(j) 


1 


S'ljj 


S9 


—S'ljj 


1 


0 


S(j) 


-S9 


-S(j) 


0 


Sip 


S9 


—S'lp 


0 



( 6 ) 



R comprises the identity plus a part equal to the (negative) cross-product 
operator [—SEx], where SE = [S'ip,69,S(j)], the vector of Euler angles ordered 
with the axes [x,y,z]. Small rotations are completely decoupled; the order of 
the small rotations does not matter. Since R~^ = R^, we have also R~^ = 
-^3x3 + SEx ; 



Xb = X — 5E X X (7) 

X = Xb + SE X Xb- (8) 

We now £x the point of interest on the body, instead of in inertial space, 
calling its location in the body frame r (radius). The differential rotations 
occur over a time step St, so that we can write the location of the point before 
and after the rotation, with respect to the first frame as follows: 



x(t) = r 

x{t + St) = R^r = r + SE x r. 



( 9 ) 



Dividing by the differential time step gives 
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Sx 5E 
It ^ It 

= u X r, 



( 10 ) 



where the rotation rate vector u ~ dE/dt because the Euler angles for this 
inhnitesimal rotation are small and decoupled. This same cross-product rela- 
tionship can be derived in the second frame as well; 



Xb{t) 

Xb{t -h St) 



Rf = f— SE X r 

r. 



such that 



( 11 ) 



Sxb SE 
It ^ It 

= U X f, 

On a rotating body whose origin point is hxed, the time rate of change of a 
constant radius vector is the cross-product of the rotation rate vector cn and 
the radius vector itself. The resultant derivative is in the moving body frame. 
In the case that the radius vector changes with respect to the body frame, we 
need an additional term: 




dxb ^ , df 

Finally, allowing the origin to move as well gives 

dxb ^ , dxo .. 

dt dt dt 

This result is often written in terms of body-referenced velocity v: 

dr 

v = uxr+ — +Vo, (15) 

where Vo is the body-referenced velocity of the origin. The total velocity of the 
particle is equal to the velocity of the reference frame origin, plus a component 
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due to rotation of this frame. The velocity equation can be generalized to any 
body-referenced vector /: 



df 

dt 



dl 

dt 



+ UJ X f. 



(16) 



2.3 Rate of Change of Euler Angles 

Only for the case of inhnitesimal Euler angles is it true that the time rate 
of change of the Euler angles equals the body-referenced rotation rate. For 
example, with the sequence [yaw, pitch, roll], the Euler yaw angle (applied hrst) 
is dehnitely not about the hnal body yaw axis; the pitch and roll rotations 
moved the axis. An important part of any simulation is the evolution of the 
Euler angles. Since the physics determine rotation rate u, we seek a mapping 
u dE/dt. 

The idea is to consider small changes in each Euler angle, and determine the 
effects on the rotation vector. The first Euler angle undergoes two additional 
rotations, the second angle one rotation, and the hnal Euler angle no additional 
rotations: 



1 


\ “ 1 


1 1 


1 1 ] 


1 1 


f S 1 


ch = R{'d)R{9)\ 




+ W) 


1 dt 1 


+ 






d<p 

t dt ) 


1 1 


i 0 J 


1 1 


[ 0 j 



1 

0 

0 



0 

cos'd 

— sin'^ 



— sin 6 
sin Ip cos 6 
cos 'd cos 9 




Taking the inverse gives 



(17) 



dE 

dt 



1 sin -d tan 6 
0 cos'd 
0 sill'd / cos 9 

T{E)u. 



cos 'd tan 9 
— sin'^ 
cos'd! cos 9 



UJ 



(18) 



Singularities exist in T at 0 = {7t/2, 37t/2}, because of the division by cos 9, and 
hence this otherwise useful equation for propagating the angular orientation of 
a body fails when the vehicle rotates about the intermediate y-axis by ninety 




2.4 Dead Reckoning 
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degrees. In applications where this is a real possibility, for example in orbiting 
satellites and robotic arms, quaternions provide a seamless mapping. For most 
ocean vessels, the singularity is acceptable, as long as it is not on the yaw axis! 

2.4 Dead Reckoning 

The measurement of heading and longitudinal speed gives rise to one of the 
oldest methods of navigation: dead reckoning. Quite simply, if the estimated 
longitudinal speed over ground is U, and the estimated heading is 0, ignoring 
the lateral velocity leads to the evolution of Cartesian coordinates: 



X = U COS0 

y = U sin0. 

Needless to say, currents and vehicle sideslip will cause this to be in error. 
Nonetheless, some of the most remarkable feats of navigation in history have 
depended on dead reckoning. 
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3 VESSEL INERTIAL DYNAMICS 

We consider the rigid body dynamics with a coordinate system affixed on the 
body. A common frame for ships, submarines, and other marine vehicles has 
the body- referenced x-axis forward, y-axis to port (left), and 2 :-axis up. This 
will be the sense of our body-referenced coordinate system here. 

3.1 Momentum of a Particle 

Since the body moves with respect to an inertial frame, dynamics expressed in 
the body-referenced frame need extra attention. First, linear momentum for 
a particle obeys the equality 



F = ^ (mv) (19) 

A rigid body consists of a large number of these small particles, which can be 
indexed. The summations we use below can be generalized to integrals quite 
easily. We have 



Fi + Ri = ^ {rriiVi ) , ( 20 ) 

where Fj is the external force acting on the particle and Ri is the net force 
exerted by all the other surrounding particles (internal forces). Since the 
collection of particles is not driven apart by the internal forces, we must have 
equal and opposite internal forces such that 

N 

Y,R. = o. ( 21 ) 

i=l 

Then summing up all the particle momentum equations gives 

^ ^ ^ d 

Y = Y W ■ (22) 

i=i i=i 

Note that the particle velocities are not independent, because the particles are 
rigidly attached. 

Now consider a body reference frame, with origin 0, in which the particle i 
resides at body-referenced radius vector r; the body translates and rotates, 
and we now consider how the momentum equation depends on this motion. 




3.2 Linear Momentum in a Moving Frame 



21 




Figure 2: Convention for the body-referenced coordinate system on a vessel: 
X is forward, y is sway to the left, and 2 : is heave upwards. Looking forward 
from the vessel bridge, roll about the x axis is positive counterclockwise, pitch 
about the y-axis is positive bow-down, and yaw about the 2 :-axis is positive 
turning left. 



3.2 Linear Momentum in a Moving Frame 

The expression for total velocity may be inserted into the summed linear mo- 
mentum equation to give 

N 

T.F. = 

i=l 



where m = and Vi = Vo+co x -Tj. Further dehning the center of gravity 

vector tq such that 



^ d 

j,{mi{vo + uj X fi)) 

i=l 

dvo d ^ 

m — — h T a; X > miTi 

dt dt ^ * 

L 2=1 



(23) 



N 

mrc = ^miri, (24) 

2=1 

we have 



dvo d ^ . 

^Fi = m— + m — {uxrG). (25) 

2 = 1 

Using the expansion for total derivative again, the complete vector equation 
in body coordinates is 
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^ , X / dvo ^ ^ did ^ 

F = ^N = m i — + u X Vo + -^ X re + uj X {u X re) V (26) 



2 = 1 



Now we list some conventions that will be used from here on: 



Vo = {u, V, w} (body-referenced velocity) 
rc = {xg, l/G, 2 :^} (body-referenced location of center of mass) 
dJ = {p, Q , (rotation vector, in body coordinates) 

F = {X,Y,Z} (external force, body coordinates). 

The last term in the previous equation simplihes using the vector triple product 
identity 



u X {u X To) = (cJ ■ rG)uJ — (iu ■ uj)fG, 

and the resulting three linear momentum equations are 



X 

Y 

Z 



dti d/Q dv 

m — + qw~rv + -^zg - -^Vg + {qVc + tzg)p - (g^ + r'^)xG 
dv dv dp 

m —+ru~pw + —xg - ~^zg + {tzg + pxG)q - {V + p^)hg 

^ ^ +pv -qu + ^vg - ^Xg + {pxG + qycY - (p^ + q^)zG 



27) 



Note that about half of the terms here are due to the mass center being in a 
different location than the reference frame origin, i.e., tg 7 ^ 0 - 



3.3 Example: Mass on a String 

Consider a mass on a string, being swung around around in a circle at speed 
U, with radius r. The centrifugal force can be computed in at least three 
different ways. The vector equation at the start is 



- f dvo ^ ^ duj _ ^ ^ ^ , 

F = m\ + u X Vo+ -- X Tg + OJ X iuj X Tg) 
\ ot dt 
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3.3.1 Moving Frame Affixed to Mass 

Affixing a reference frame on the mass, with the local x oriented forward and 
y inward towards the circle center, gives 



such that 



u 

re 

dvo 

dt 

ddo 



{0,0,U/rf 

{ 0 , 0 , 0 }^ 

{ 0 , 0 , 0 }^ 

{ 0 , 0 , 0 }^, 



F = mu X Vo = m{0, U‘^/r, 0}^. 

The force of the string pulls in on the mass to create the circular motion. 

3.3.2 Rotating Frame Attached to Pivot Point 

Affixing the moving reference frame to the pivot point of the string, with the 
same orientation as above but allowing it to rotate with the string, we have 



u 

re 

dvp 

dt 

du 

giving the same result: 



{ 0 , 0 , 0 }^ 

{0,r,0}^ 

{ 0 , 0 , 0 }^ 

{ 0 , 0 , 0 }^, 



F = mu X {u X tg) = m{0, U‘^/r, 0}^. 
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3.3.3 Stationary Frame 

A frame fixed in inertial space, and momentarily coincident with the frame 
on the mass (3.3.1), can also be used for the calculation. In this case, as the 
string travels through a small arc Sip, vector subtraction gives 

Sv = {0, t/ sin (5-^,0}'^ ~ {0, 1/(5-^, 0}^. 

Since ip = U/r, it follows easily that in the hxed frame dv/dt = {0, /r, 0}^, 

as before. 



3.4 Angular Momentum 

For angular momentum, the summed particle equation is 
^ ^ ^ ^ d 

Y^{Mi + fi X Fj) = X —iniiVi), (28) 

i=l i=l 

where Mj is an external moment on the particle i. Similar to the case for linear 
momentum, summed internal moments cancel. We have 



N ^ ^ ^ 

Y^{Mi + TiX Fi) = rmn X 

i=l i=l 

N 

X (dj X (C X fj)). 

i=l 

The summation in the hrst term of the right-hand side is recognized simply as 
mrc, and the hrst term becomes 



dvo 

ht 



+ id X Vo 




xfi] + 



mrc X 



dvo 

ht 



d X Vo . 



The second term expands as (using the triple product) 



(29) 



N 

Y ^ 

i=l 




N 

i=l 






f + zf)p- {yiq + Zir)xi) j 

I EZi + zZq - Zip + Zif)yi) \ 
y EZi iZf + yl)r - {Xip + yiq)zi) J 



(30) 
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Employing the definitions of moments of inertia, 



I 



'■yy 



•■xy 



•■yz 



^xx ^xy ^xz 

lyx lyy lyz 

L Izx Izy Izz J 
N 



miiVi + z^) 

i=l 

N 

E + z^) 

i=l 

N 

E + Vi) 

i=l 

N 

lyx ^ ' 'IXliXigi 

i=l 

N 

hx = ~Y1 

i=l 

N 

Izy ^ ' xriiyiZi.j 



i=l 



(inertia matrix) 



(cross-inertia) 



the second term of the angular momentum right-hand side collapses neatly 
into Idu/dt. The third term can be worked out along the same lines, but 
offers no similar condensation: 



N 

mifi X ((c5 ■ fi)Q - {Q ■ Q)fi) 

i=l 



^ ruifi X c5(ci; ■ fj) (31) 

i=l 

{ - Ziq){xip + Viq + Ziv) 1 

< Y.i=imi{zip - Xir){xip + yiq + Ziv) > 

I Y.f=imi{xiq - yip){xip + yiq + Zir) J 

f Iyz{q^ - r^) -f I^^pq - I^ypr j 

< Ixzir"^ -P^) + Ixyrq - lyzPq > + 

[ Ixy{p^ - q^) + lyzpr - Ixzqr ) 

( {Izz ~ Iyy)'I^Q I 

\ {^Ixx Izz'l^P / • 

I ^^yy ~ Ixx)qp ) 
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Letting M = {K, M, N} be the total moment acting on the body, i.e., the left 
side of Equation 28, the complete moment equations are 

K IxxP IxyQ A Ixzl A (2^) 

{hz - Iyy)rq A Iyz{q^ ~ + RzPq ~ hyPr + 

m [|/g(w + pv — qu) — zg{v + ru — pw)] 

^yxP ^yyQ 

{Ixx - Izz)pr A Ixz{r^ - P^) A Ryqr - ly^qp A 
m [zciii A qw — rv) — xc{w A pv — qu)] 

N IzxP A Izyq A IzzV A 

i-^yy Zx)pq A Ixyi^P q ) A lyzpv Ixzqv ~l“ 
m [xg{v a ru — pw) — yciu A qw — rv)] . 

3.5 Example: Spinning Book 

Consider a homogeneous rectangular block with I^x < lyy < Rz and all off- 
diagonal moments of inertia are zero. The linearized angular momentum equa- 
tions, with no external forces or moments, are 

Ixx T i^Izz Iyy')^q 0 

lyy T (Lxx Izz)pV 0 

dr 

Izz~^ A \Iyy Ixx)qP 0- 

We consider in turn the stability of rotations about each of the main axes, 
with constant angular rate The interesting result is that rotations about 
the X and ^ axes are unstable, while rotation about the y axis is not. 

3.5.1 x-axis 

In the case of the x-axis, p = Q A Sp, q = Sq, and r = Sr, where the S prehx 
indicates a small value compared to 12. The hrst equation above is uncoupled 
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from the others, and indicates no change in 6p, since the small term 6q6r can 
be ignored. Differentiate the second equation to obtain 



. d^6q 
yy-^ 



B i^Ixx 






d6r 



0 



Substitution of this result into the third equation yields 



d‘^Sq 



lyyhz ^^2 



i^lxx kzz){^lxx - Iyy)kf5q 



0 . 



A simpler expression is 5q + a5q = 0, which has response 5q{t) = 
when (5g(0) = 0. For spin about the x-axis, both coefficients of the differential 
equation are positive, and hence a > 0. The imaginary exponent indicates 
that the solution is of the form 5q{t) = 5q{i5)cos^/oit, that is, it oscillates but 
does not grow. Since the perturbation 5r is coupled, it too oscillates. 



3.5.2 y-axis 

Now suppose q = kl + 5q\ differentiate the first equation and substitute into 
the third equation to obtain 



I ZZ I XX 



d‘^Sp 

~w 



^^yy ^x:x^(.^yy ^zz^^ 



0 . 



Here the second coefficient has negative sign, and therefore a < 0. The ex- 
ponent is real now, and the solution grows without bound, following Sp(t) = 



3.5.3 2 :-axis 

Finally, let r = D -|- 5r: differentiate the first equation and substitute into the 
second equation to obtain 



d'^Sp 



lyyixx ^^2 



T i^^xx Izz^i^Iyy ^zz')^ ^P 



0 . 



The coefficients are positive, so bounded oscillations occur. 




3 VESSEL INERTIAL DYNAMICS 



3.6 Parallel Axis Theorem 

Often, the mass center of an body is at a different location than a more con- 
venient measurement point, the geometric center of a vessel for example. The 
parallel axis theorem allows one to translate the mass moments of inertia refer- 
enced to the mass center into another frame with parallel orientation, and vice 
versa. Sometimes a translation of coordinates to the mass center will make the 
cross-inertial terms Ry, lyz, Rz small enough that they can be ignored; in this 
case r G = 0 also, so that the equations of motion are signihcantly reduced, as 
in the spinning book example. 

The formulas are: 

Ixx = Lx + m{5y^ + 5zL (33) 

lyy = lyy Y Tfl{5X “|- 5Z ) 

Lz = Lz + rn{5x^ + 5yL 
Lz = Lz - m6y6z 
Lz = Lz — mSx6z 
Ly Ly TTldxdy ^ 

where I represents an MMOI in the axes of the mass center, and 6x, for ex- 
ample, is the translation of the x-axis to the new frame. Note that translation 
of MMOI using the parallel axis theorem must be either to or from a frame 
resting exactly at the center of gravity. 

3.7 Basis for Simulation 

Except for external forces and moments F and M, we now have the neces- 
sary terms for writing a full nonlinear simulation of a rigid body, in body 
coordinates. There are twelve states, comprising the following components: 

• Vo, the vector of body-referenced velocities. 

• uj, body rotation rate vector. 

• X, location of the body origin, in inertial space. 

• E, Euler angle vector. 
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The derivatives of body-referenced velocity and rotation rate come from Equa- 
tions 27 and 32, with some coupling which generally requires a 6 x 6 matrix 
inverse. The Cartesian position propagates according to 

X = {E)vo, (34) 

while the Euler angles follow: 

E = V{E)u. (35) 
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4 HYDRODYNAMICS: INTRODUCTION 

The forces and moments on a vessel are complicated functions of many fac- 
tors, including water density, viscosity, surface tension, pressure, vapor pres- 
sure, and motions of the body. The most important factors for large ocean 
vehicles are density and motion, and we can make simplihcations to param- 
eterize the most prominent relationships. This section pertains to the use of 
hydrodynamic coefficients for predicting hydrodynamic response. 

4.1 Taylor Series and Hydrodynamic Coefficients 

Recall the Taylor expansion of a function: 

f(x) = f{xo) + - ^o) + + ■■■■ (36) 

We introduce the notation 

dfjxo) 

HT 

1 

2 ! dx'^ ’ 

and so on, so that a two-variable Taylor expansion would have the form 

f{x,y) = f{xo,yo)+ (37) 

fx(x-Xo) + fyiy-yo) + 

fxx{x Xq) T fyyiV Vo) T fxy{x Xo){y |/o) T 
fxxx{x - Xof H . (38) 

Note that all of the factorials are included in the coefficients. This notation 
covers some instances where the formal Taylor series is meaningless, but the 
notation is still clear. As one example, fluid drag is often written as 

F = ^pCdAu\u\ = Fu\u\u\u\. 



fx = 

fxX = 




4.2 Surface Vessel Linear Model 
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4.2 Surface Vessel Linear Model 

We now discuss some of the hydrodynamic parameters which govern a ship 
maneuvering in the horizontal plane. The body x-axis is forward and the y- 
axis is to port, so positive r has the boat turning left. We will consider motions 
only in the horizontal plane, which means 6 = ijj = p = q = w = 0. Since the 
vessel is symmetric about the x — z plane, yc = 0; zq is inconsequential. We 
then have at the outset 



X = 
Y = 
N = 

Letting u = U + u, where U 
set is 



m 



' du 

A 

' dv 



rv — xcr 



(39) 



dr\ 

dr f dv \ 



» u, and eliminating higher-order terms, this 



X 

Y 

N 



m- 



( TT 

m \ — + rU + a 
\ dt 

^ dr ( c 



(9r\ 

+ rU 



(40) 



A number of coefficients can be discounted. First, in a homogeneous sea, with 
no current, wave, or wind effects, {X^, Xy, Xfi,,Yx,Yy,Y^, N^, Ny, N^} are all 
zero. We assume that no hydrodynamic forces depend on the position of the 
vessel.^ Second, consider X„: since this longitudinal force would have the same 
sign regardless of the sign of v (because of side-to-side hull symmetry), it must 
have zero slope with v at the origin. Thus X.^ = 0. The same argument shows 
that {Xr, Xi,, Xj.,Yu,Yu, Nu, Nu} = 0. Finally, since fluid particle acceleration 
relates linearly with pressure or force, we do not consider nonlinear acceleration 

^Note that the linearized heave/pitch dynamics of a submarine do depend on the pitch 
angle; this topic will be discussed later. 
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terms, or higher time derivatives. It should be noted that some nonlinear terms 
related to those we have eliminated above are not zero. For instance, Yuu = 0 
because of hull symmetry, but in general Xy^ = 0 only if the vessel is bow-stern 
symmetric. 

We have so far, considering only the linear hydrodynamic terms. 



[m 



(m - Xy)ii 


= XyUPX' 


(41) 


- Yy)v + {mxG — Yy)r 


= YyV + {Yr-mU)r + Y' 


(42) 


- Ny)v + - Ny)r 


= NyV — {Ny — mxGU)r + N' . 


(43) 



The right side here carries also the imposed forces from a thruster(s) and 
rudder(s) {X' ,Y' , N'}. Note that the surge equation is decoupled from the 
sway and yaw, but that sway and yaw themselves are coupled, and there- 
fore are of immediate interest. With the state vector s = {u,r} and exter- 
nal force/moment vector F = {Y',N'}, a state-space representation of the 
sway/yaw system is 



m — Yy 


mxG — Yy 


ds 


■ Yy 


Yy-mU 


mxG — Ny 


IzZ-Ny 


dt 


_ Ny 


Ny — uixgU 



Ms = Ps + F 

k = M-^Ps + M~^F 



(44) 



s = As + BF. 



(45) 



The matrix M is a mass or inertia matrix, which is always invertible. The 
last form of the equation is a standard one wherein A represents the internal 
dynamics of the system, and B is a gain matrix for the control and disturbance 
inputs. 



4.3 Stability of the Sway /Yaw System 

Consider the homogeneous system s = As: 



S\ 



-|- A12S2 
^21'Sl + ^22"S2- 
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We can rewrite the second equation as 

~ ^22^ ^21'Sl 

and substitute into the first equation to give 

Si + (—^11 — ^22)Sl + {A 11 A 22 — Ai2^2i)si = 0. (47) 

Note that these operations are allowed because the derivative operator is linear; 
in the language of the Laplace transform, we would simply use s. A necessary 
and sufficient condition for stability of this ODE system is that each coefficient 
must be greater than zero: 




(46) 



—All — A22 > 0 

A11A22 A12A21 > 0 

The components of A for the sway/yaw problem are 



(48) 



{Izz - Nr)Y^ + (Yf. - mxG)Ny 

(m - W)(A^ - Nr) - {mxG - Yf){mxG - N^,) 

— {Izz — Nr) {mil — Yr) — {Yr — mXG){TnXGU — Nr) 

{m - Yr){Izz - Nr) - {mXG - Yr){mXG - Ni,) 

{Nr - mXG)Yr + {m- Yj,)Nr 

(m - Yi,){Izz - Nf) - {mxG - Yr){mxG - Ni,) 

— {Ni, — mxG){mU — Yr) — {m — Yi,){mxGU — W) 
(m - Yi,){Izz - Nr) - {mxG - Yr){mxG - Ni,) 



The denominators are identical, and can be simplihed. First, let xg — 0; 
valid for many vessels with the origin is at the geometric center. If the origin 
is at the center of mass, xg = 0. Next, if the vessel is reasonably balanced 
with regard to forward and aft areas with respect to the origin, the terms 
{Ni,, Yr, Nr, Yr} take very small values in comparison with the others. To wit, 
the added mass term —17 is of the order of the vessel’s material mass m, 
and similarly W ~ —Izz- Both 17 and A7 take large negative values. Linear 
drag and rotational drag are signihcant also; these are the terms 17 and Nr, 
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both large and negative. The denominator for A’s components reduces to 
(m - Yi,){R^ - Nr), and 



^11 

A22 



Yr 

m—Yf, 

Nr 

hz-Nr 



< 0 
< 0 . 



Hence the hrst condition for stability is met: —Hu — A22 > 0. For the second 
condition, since the denominators of the Aij are identical, we have only to look 
at the numerators. For stability, we require 



- Nr)Yr{m - Yr)Nr (50) 

- [NrYr + (m - Y^Nr] [-(/,, - Nr){mU - Yr) + YrNr] > 0. 

The hrst term is the product of two large negative and two large positive 
numbers. The second part of the second term contains mU, which has a large 
positive value, generally making stability critical on the (usually negative) N^. 
When only the largest terms are considered for a vessel, a simpler form is 
common: 



C = YrNr + NrimU - Yr) > 0. (51) 

C is called the vessels stability parameter. The terms of C compete, and 
yaw/sway stability depends closely on the magnitude and sign of W- Adding 
more surface area aft drives Ny more positive, increasing stability as ex- 
pected. Stability can also be improved by moving the center of gravity forward. 
Nonzero xg shows up as follows: 

C = Yy{Nr-mxGU) + Ny{mU -Yr)>0. (52) 

Since W and Yy are both negative, positive xg increases the (positive) inhuence 
of C’s hrst term. 

4.4 Basic Rudder Action in the Sway/ Yaw Model 

Rudders are devices which develop large lift forces due to an angle of attack 
with respect to the oncoming huid. As in our discussion of lift on the body, 
the form is as follows: L = ^pAU^Ci(a), where a is the angle of attack. The 




4.4 Basic Rudder Action in the Sway/Yaw Model 



35 




Figure 3: Convention for positive rudder angle in the vessel reference system. 



lift coefficient Ci is normally linear with a near a = 0, but the rudder stalls 
when the angle of attack reaches a critical value, and thereafter develops much 
less lift. We will assume that a is small enough that the linear relationship 
applies: 



Q{a) 



dQ 

da 



a. 

q :=0 



(53) 



Since the rudder develops force (and a small moment) far away from the body 
origin, say a distance I aft, the moment equation is quite simple. We have 









Q :=0 



IU\ 



q ;=0 



(54) 

(55) 



Note the difference between the rudder angle expressed in the body frame, 6, 
and the total angle of attack a. Angle of attack is influenced by S, as well 
as v/U and Ir. Thus, in tank testing with v = 0, S = a and Ns = Na, 
etc., but in real conditions, other hydrodynamic derivatives are augmented to 
capture the necessary effects, for example W and W- Generally speaking, the 
hydrodynamic characteristics of the vessel depend strongly on the rudder, even 
when 5 = 0. In this case the rudder still opposes yaw and sway perturbations 
and acts to stabilize the vessel. 

A positive rudder deflection (defined to have the same sense as the yaw angle) 
causes a negative yaw perturbation, and a very small positive sway perturba- 
tion. 
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4.4.1 Adding Yaw Damping throngh Feedback 

The stability coefficient C resulting from the addition of a control law S = krV, 
where > 0 is a feedback gain, is 

C = Yy{Nr — rrixcU + kyNs) + Ny{mU — Yy — KYs). (56) 

Ys is small positive, but Ns is large and negative. Hence C becomes more 
positive, since Yy is negative. 

Control system limitations and the stalling of rudders make obvious the fact 
that even a very large control gain ky cannot completely solve stability prob- 
lems of a poorly-designed vessel with an inadequate rudder. On the other 
hand, a vessel which is overly stable (C* >> 0 with no rudder action) is unma- 
neuverable. A properly-balanced vessel just achieves stability with zero rudder 
action, so that a reasonable amount of control will provide good maneuvering 
capabilities. 



4.4.2 Heading Control in the Sway/Yaw Model 

Considering just the yaw equation of motion, i.e., n = 0, with a rudder, we 
have 



(4, - Ny)^ + {mxGU - Ny)^ = NsS. (57) 

Employing the control law S = k^pcj), the system equation becomes a homoge- 
neous, second-order ODE: 



- Ny)^ -f {mxcU - Ny)^ - Nsk^4> = 0. (58) 

Since all coefficients are positive (recall < 0), the equation gives a stable 9 
response, settling under second-order dynamics to 9{oo) = 0. The control law 
S = k^{4) — (pdesired) + kyV is the basis for heading autopilots, which are used 
to track (pdesired- This use of an error signal to drive an actuator is in fact the 
essence of feedback control. In this case, we require sensors to obtain r and 0, 
a controller to calculate S, and an actuator to implement the corrective action. 
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4.5 Response of the Vessel to Step Rudder Input 

4.5.1 Phase 1: Accelerations Dominate 

When the rudder first moves, acceleration terms dominate, since the velocities 
are zero. The equation looks like this: 



m — Yy mxG — Yf. 
mxG - Ni, - Nr 

Since W and mxG are comparatively small in the first row, we have 

t.(0) = (60) 

m — Yj, 

and the vessel moves to the left, the positive n-direction. The initial yaw is in 
the negative r-direction, since Ns < 0: 




f(0) = 



Ns6 

hz-NR 



(61) 



The hrst phase is followed by a period (Phase 2), in which many terms are 
competing and contributing to the transient response. 



4.5.2 Phase 3: Steady State 

When the transients have decayed, the vessel is in a steady turning condition, 
and the accelerations are zero. The system equations reduce to 

n 1 _ ^ f {tuxgU - Nr)Ys + (W - mU)Ns 
r j~ C\ NrYs - YrNs 

Note that the denominator is the stability parameter. The steady turning rate 
is thus approximated by 




r = 



Y,Ns 



C 



-s. 



(63) 



With C > 0, the steady-state yaw rate is negative. If the vessel is unstable 
(C* < 0), it turns in the opposite direction than expected. This turning rate 
equation can also be used to estimate turning radius R: 



R = 



U 



uc 

-YrNs5' 



r 



(64) 
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The radius goes up directly with C, indicating that too stable a ship has poor 
turning performance. We see also that increasing the rudder area increases 
Ns, decreasing R as desired. Increasing the deflection 6 to reduce R works 
only to the point of stalling. 

4.6 Summary of the Linear Maneuvering Model 

We conclude our discussion of the yaw/sway model by noting that 

1. The linearized sway/yaw dynamics of a surface vessel are strongly cou- 
pled, and they are independent of the longitudinal dynamics. 

2. The design parameter C should be slightly greater than zero for easy 
turning, and “ hands-off” stability. The case C < 0 should only be 
considered under active feedback control. 

3. The analysis is valid only up to small angles of attack and turning rates. 
Very tight maneuvering requires the nonlinear inertial components and 
hydrodynamic terms. Among other effects, the nonlinear equations cou- 
ple surge to the other motions, and the actual vessel loses forward speed 
during maneuvering. 

4.7 Stability in the Vertical Plane 

Stability in the horizontal plane changes very little as a function of speed, 
because drag and lift effects generally scale with U"^. This fact is not true 
in the vertical plane, for which the dimensional weight/buoyancy forces and 
moments are invariant with speed. For example, consider the case of heave 
and pitch, with xg = 0 and no actuation: 



Z^w ZyjW Zqq -\- Zqq -\- i^B — W) (65) 

MyjW -|- MyjW + Mqq + Mqq — Bis sin 6. (66) 

The last term in each equation is a hydrostatic effect induced by opposing net 
buoyancy B and weight W . U denotes the vertical separation of the center 
of gravity and the center of buoyancy, creating the so-called righting moment 
which nearly all underwater vehicles possess. Because buoyancy effects do not 



m 



Uq\ = 
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change with speed, the dynamic properties and hence stability of the vehicle 
may change with speed. 
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5 SIMILITUDE 

5.1 Use of Nondimens ional Groups 

For a consistent description of physical processes, we require that all terms 
in an equation must have the same units. On the basis of physical laws, 
some quantities are dependent on other, independent quantities. We form 
nondimensional groups out of the dimensional ones in this section, and apply 
the technique to maneuvering. 

The Buckingham 7r-theorem provides a basis for all nondimensionalization. 
Let a quantity Qn be given as a function of a set of n — 1 other quantities: 

Qn = fQ{Ql,Q2,- ■ ■ ,Qn-l)- (67) 

There are n variables here, but suppose also that there are only k independent 
ones; k is equivalent to the number of physical unit types encountered. The 
theorem asserts that there are n — k dimensionless groups tt* that can be 
formed, and the functional equivalence is reduced to 

'^n—k ^2) ' ' ' i fc— l)- (68) 

Example. Suppose we have a block of mass m resting on a frictionless hor- 
izontal surface. At time zero, a steady force of magnitude F is applied. We 
want to know X{T), the distance that the block has moved as of time T. The 
dimensional function is X(T) = fQ^m, F,T), so n = 4. The (MKS) units are 



[^(■)] 


= m 


[m] 


= kg 


[F] 


= kgm/s 


|T] 


= s, 



and therefore k = 3. There is just one nondimensional group in this re- 
lationship; 7Ti assumes only a constant (but unknown) value. Simple term- 
cancellation gives 7Ti = X{T)m/FT‘^, not far at all from the known result that 
X{T) = FT^/2m\ 

Example. Consider the flow rate Q of water from an open bucket of height 
h, through a drain nozzle of diameter d. We have 
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Q = fQ{h,d,p,^,g), 

where the water density is p, and its absolute viscosity p; g is the acceleration 
due to gravity. No other parameters affect the flow rate. We have n = 6, and 
the (MKS) units of these quantities are: 



[Q] = rn^/s 
[h] = m 
[d\ = m 
[p] = kg/m^ 

[fj] = kg /ms 

[g\ = m/ 

There are only three units that appear: [length, time, mass], and thus k = 3. 
Hence, only three non-dimensional groups exist, and one is a unique func- 
tion of the other two. To arrive at a set of valid groups, we must create 
three nondimensional quantities, making sure that each of the original (di- 
mensional) quantities is represented. Intuition and additional manipulations 
come in handy, as we now show. 

Three plausible hrst groups are: tti = pQ/dp, 712 = dp^/gh/p,, and tts = h/d. 
Note that all six quantities appear at least once. Since h and d have the same 
units, they could easily change places in the hrst two groups. However, tti is 
recognized as a Reynolds number pertaining to the orihce how. 7T2 is more 
awkward, but products and fractions of groups are themselves valid groups, 
and we may construct 7T4 = = Q/d^y/gh to nondimensionalize Q with 

a pressure velocity, and then tts = 7ri/7T4 = pd^/gh/ p, to establish an orihce 
Reynolds number independent of Q. We hnally have the useful result 



^4 /7 t (^57^2) ^ 

Q ^ f ( h\ 

d‘^y/gh \ p ' dj 

The uncertainty about where to use h and d, and the questionable importance 
of h/d as a. group are remnants of the theorem. Intuition is that h/d is im- 
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material, and the other two terms have a nice physical meaning, e.g., tts is a 
Reynolds number. 

The power of the 7r-theorem is primarily in reducing the number of parameters 
which must be considered independently to characterize a process. In the flow 
example, the theorem reduced the number of independent parameters from 
hve to two, with no constraints about the actual physics taking place. 



5.2 Common Groups in Marine Engineering 

One frequently encounters the following groups in fluid mechanics and marine 
engineering: 



1. Fronde number: 



Fr = 



U 



v^’ 



(69) 



where U is the speed of the vessel, g is the acceleration due to gravity, 
and L is the waterline length of the vessel. The Fronde number appears 
in problems involving pressure boundary conditions, such as in waves on 
the ocean surface. Roughly speaking, it relates the vessel speed U to 
water wave speeds of wavelength L; the phase speed of a surface wave is 
V = \g/2Ti, where A is the wavelength. 



2. Cavitation number: 



p 

5= ^ 



F, 



yu^ ’ 



(70) 



where Pqo represents the ambient total pressure, the vapor pressure 
of the fluid, and U the propeller inlet velocity. A low cavitation number 
means that the Bernoulli pressure loss across the lifting surface will cause 
the fluid to vaporize, causing bubbles, degradation of performance, and 
possible deterioration of the material. 



3. Reynolds number: 



Re 



Ul 

i^/p’ 



(71) 




5.2 



Common Groups in Marine Engineering 



43 



where U is velocity, p is absolute viscosity, and p is density. Since Re 
appears in many applications, I represents one of many length scales. 
Reynolds number is a ratio of fluid inertial pressures to viscous pres- 
sures: When Re is high, viscous effects are negligible. Re can be used to 
characterize pipe flow, bluff body wakes, and flow across a plate, among 
others. 

4. Weber number: 



W = 




( 72 ) 



where a is the surface tension of a fluid. Given that [a] = N/m (MKS), 
p[/^ normalizes pressure, and I normalizes length. The Weber number 
indicates the importance of surface tension. 



To appreciate the origins of these terms from a fluid particle’s point of view, 
consider a box having side lengths [dx, dy, dz]. Various forces on the box scale 
as 



(inertia) Fi 

(gravity) Fg 
(pressure) Fp 

(shear) Fg 
(surface tension) F^ 



dv dv 

= p—dxdydz -\- pv—dxdydz ~ pU“^P 
at ox 

= pgdxdydz ~ pgl^ 

= Pdxdy ~ Pf 

= p^dxdy ~ p,Ul 
dz 

= adx ~ al. 



Thus the groups listed above can be written as 



Fr = 



Fi 






F, 


gi 


Fr, 


P 






Re = ^ 


pUl 


Fs 


y 


W = — 


pUH 


Frr 


a 
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When testing models, it is imperative to maintain as many of the nondimen- 
sional groups as possible of the full-scale system. This holds for the geometry 
of the body and the kinematics of the flow, the surface roughness, and the 
all of the relevant groups governing fluid dynamics. Consider the example of 
nozzle flow from a bucket. Suppose that we conduct a model test in which 
Re is abnormally large, i.e., the viscous effects are negligible. Under inviscid 
conditions, the flow rate is Q = 7rd‘^y/2gh/4. This rate cannot be achieved for 
lower-i?e conditions because of fluid drag in the orifice, however. 

In a vessel, we write the functional relationship for drag as a starting point: 



Cr 



D 

\pAU^ 



U{Re,Fr,W). 



First, since I is large, W is very large, and hence surface tension plays no role. 
Next, we look at Re = Ulju and Fr = U/y/^, both of which are important 
for surface vessels. Suppose that Uhip = ^Imodeh so that usually A >> 1; 
additionally, we set gmodei = 9ship-, i-e., the model and the true vessel operate 
in the same gravity held. 

Fronde number similitude requires Umodei = Ughip/ V^- Then Reynolds number 
scaling implies directly i^modei = Unfortunately, few huids with this 

property are workable in a large testing tank. As a result, accurate scaling of 
Re for large vessels to model scale is quite difficult. 

For surface vessels, and submarines near the surface, it is a routine procedure 
to employ turbulence stimulators to achieve how that would normally occur 
with ship-scale Re. Above a critical value Re ~ 500,000, Cf is not sensitive 
to Re. With this achieved, one then tries to match Fr closely. 



5.3 Similitude in Maneuvering 

The linear equations of motion for the horizontal yaw/sway problem are: 



(m — Yi,)v — Y^v -f {mU — Yr)r + {mxc — Yf.)r = Y 
{Rz — Nry + {mxcU — Nr)r + {rrixc — NRi) — N^v = N. 



These equations can be nondimensionalized in a standard way, by using the 
quantities [U,L,p]: these three values provide the necessary units of length. 
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time, and mass, and furthermore are readily accessible to the user. First, we 
create nondimensional states, denoted with a prime symbol; 




We follow a similar procedure for the constant terms as follows, including a 
factor of 1/2 with p, for consistency with our previous expressions: 



x'g 

U' 

Y'- 

Y' 

Y' 

Y' 

Y' 

N' 

j. y y 

N' 

V 



\pL^ 

^zz 

w 

L 

ipL3 

w 

IpUL^ 

Yr 

\pL^ 

Yr 

\pUL^ 

Y 

JpLf^ 

Nr 

\pL^ 

Nr 

\pUL^ 



(74) 
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M 



Ni 



N' 



Nr 

w 

Nr 

\pUL^ 

N 

JpUH?' 



Note that every force has been normalized with and every moment 

with ^pU'^L^] time has been also nondimensionalized with L/U. Thus we 
arrive at a completely equivalent set of nondimensional system equations, 



(m' -Y')E -ry + (m'U' -Y;y + (m'xG-Yiy = Y' (75) 
(/', - N'y + {m'xcU' - N'r'y + - N'y - N'y = N' . 

Since fluid forces and moments generally scale with U"^, the nondimensionalized 
description holds for a range of velocities. 

5.4 Roll Equation Similitude 

Certain nondimensional coefficients may arise which depend explicitly on U, 
and therefore require special attention. Let us carry out a similar normaliza- 
tion of the simplihed roll equation 

(4x - Kp)p - Kpp - Ky = K. (76) 

For a surface vessel, the roll moment is based on metacentric stability, and 
has the form = —pg'V{GM), where V is the displaced fluid volume of the 
vessel, and GM is the metacentric height. The nondimensional terms are 



IL 



K'p 



K 




\pL^ 

Kp 

\pL^ 

Kp 

\pUL^ 

\pUM^ 



( 77 ) 
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P = 



P = 



P 



7/2 

L 

Z7^’ 



leading to the equivalent system 

(/L - - A>' - V'(GM> = K'. (78) 

Note that the roll angle 0 was not nondimensionalized. The Fronde number 
has a very strong influence on roll stability, since it appears explicitly in the 
nondimensional righting moment term, and also has a strong influence on K^. 
In the case of a submarine, the righting moment has the form = —Bh, 
where B is the buoyant force, and h is the righting arm. The nondimensional 
coefficient becomes 



^ \pU^L^' 

again depends strongly on U, since B and h are fixed; this needs to be 
maintained in model tests. 
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6 CAPTIVE MEASUREMENTS 

Before making the decision to measure hydrodynamic derivatives, a prelimi- 
nary search of the literature may turn up useful estimates. For example, test 
results for many hull-forms have already been published, and the basic lift- 
ing surface models are not difficult. The available computational approaches 
should be considered as well; these are very good for predicting added mass 
in particular. Finally, modern sensors and computer control systems make 
possible the estimation of certain coefficients based on open-water tests of a 
model or a full-scale design. 

In model tests, the Froude number Fr = which scales the influence 

of surface waves, must be maintained between model and full-scale surface 
vessels. Reynolds number Re = — , which scales the effect of viscosity, need 
not be matched as long as the scale model attains turbulent flow (supercritical 
Re). One can use turbulence stimulators near the bow if necessary. Since the 
control surface(s) and propeller(s) affect the coefficients, they should both be 
implemented in model testing. 



6.1 Towtank 

In a towtank, tow the vehicle at different angles of attack, measuring sway 
force and yaw moment. The slope of the curve at zero angle determines 
and iV„ respectively; higher-order terms can be generated if the points deviate 
from a straight line. Rudder derivatives can be computed also by towing with 
various rudder angles. 



6.2 Rotating Arm Device 

On a rotating arm device, the vessel is hxed on an arm of length R, rotating at 
constant rate r: the vessel forward speed isU = rR. The idea is to measure the 
crossbody force and yaw moment as a function of r, giving the coefficients Yr 
and Nr- Note that the lateral force also contains the component {m — Yi,)r‘^R. 
The coefficients Yy and Ny can also be obtained by running with a hxed angle 
of attack. Finally, the measurement is made over one rotation only, so that 
the vessel does not re-enter its own wake. 




6.3 Planar-Motion Mechanism 
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6.3 Planar-Motion Mechanism 

With a planar motion mechanism, the vessel is towed at constant forward 
speed U, bnt is held by two posts, one forward and one aft, which can each 
impose independent sway motions, therefore producing variable yaw. The 
model moves in pure sway if ya(t) = ytit), in a pure yaw motion about the 
mid-length point if ya{t) = —yb{t), or in a combination sway and yaw motion. 
The connection points are a distance I forward and aft from the vessel origin. 
Usually a sinusoidal motion is imposed: 



ya(t) = a cos cut (79) 

yb(t) = b cos{ut + 'ip) , 

and the transverse forces on the posts are measured and approximated as 



Ya(t) = FaCOs{ut + 9a) (80) 

Yb(t) = Fb cos{ut + 9b) . 

If linearity holds, then 



(m — Yi,)v — Wn -|- {mil — Yr)r + {mxc — Yr)r 
{hz - Nr)r + {mxcU - Nr)r + {mxc - Ni,)v - N.^v 



Ya + Yb (81) 
{Yb - Ya)l. 



We have n = (?/« + w)/2 and r = {jjb — ya)/‘^h When a = b, these become 



V 

V 

r 

r 



aoj 

T 



(sina;t(l -f cos'^) -f cos cut sin ■^) 



~Y 



aoj 



■ (cosci;t(l -1- cosip) — sinutsimp) 
(sina;t(cos'^ ~ 1) + cosent sin ■^) 



acn^ 



(coscnt(cos'^ — 1) 



sinent sin'^) . 



(82) 



Equating the sine terms and then the cosine terms, we obtain four independent 
equations: 
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(m - n) 



au 

'T" 

V, 



(1 + COS'lp) — 



( , 

( — 2^ ) sin^ + 

( aoo\ 

[mU - Yr) sin'?/’ + 



{mxG - Yr) 



(m - Yr) 



Y, 



aoj 



{COS'lp — 1 ) 



(— siii'^) — 



(“t) + + 

{mU - Yr) {COS'lp - 1) + 



{mxG - Yr) 



aoj 






{Izz - Nr) 



aoj 



{COS'lp — 1) + 



{tUXgU — Nr) sill'^ + 

'' au‘^\ 

— I (1 +COS?/’) - 

( au\ . , 

Nvi — y) 



{mxG - Ni,) 



{h^ - Nr) (- sin?/’) + 

{ttixgU - Nr) {COS'lp - 1) + 

{mxG - Ni,) {- sin 'ip) - 

Nv (1 + COS'lp) 



(83) 



Fa COS 6a + Fb COS 6b 



(84) 



-Fa sin 6 a - Fb sin 6b 



(85) 



l{Fb cos 6b - Fa cos 6a) 



( 86 ) 



1{-Fbsin 6b + Fasin6a) 
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In this set of four equations, we know from the imposed motion the values 
[U, tjj, a, oj]. From the experiment, we obtain [F^, F;,, 6a, and from the rigid- 
body model we have [m, Izz,xg\- It turns out that the two cases of = 0 
(pure sway motion) and tjj = 180° (pure yaw motion) yield a total of eight 
independent equations, exactly what is required to hnd the eight coefficients 
\Yi,,Y.a,Yr,yr, Ni,, Ny., Nr]. Remarkably, we can write the eight solutions 
directly: For = 0, 



-n(-^)p) . 

(mxa - Ni) (2) = 

-/V.(-f)(2) = 

to be solved respectively for [Yy,Yr,Ni,, 



-- Fa COS 6a + Fb COS 6b 
-- -FaSin6a - Fbsin6b 
-- l{FbCos6b- Fa cos 6a) 

-- 1{-Fbsin6b + FaSin6a), 
iV„]. For = 180°, we have 



(87) 



[mxc - Yf) (-2) 

{mU - y;) (-K) (-2) 

{mxcU -Nr) (- 2 ) 



Fa COS 6a + Fb COS 6b (88) 

-Fa sin 6a - Fb sin 6b 
l{Fb cos 6b - Fa cos 6a) 

1{-Fbsm6b + Fa sin 6a), 



to be solved for [F^, Yr, Nf., Nr]. Thus, the eight linear coefficients for a surface 
vessel maneuvering, for a given speed, can be deduced from two tests with 
a planar motion mechanism. We note that the nonlinear terms will play a 
signihcant role if the motions are too large, and that some curve htting will be 
needed in any event. The PMM can be driven with more complex trajectories 
which will target specihc nonlinear terms. 
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7 STANDARD MANEUVERING TESTS 

This section describes some of the typical maneuvering tests which are per- 
formed on full-scale vessels, to assess stability and performance. 

7.1 Dieudonne Spiral 

1. Achieve steady speed and direction for one minute. No changes in speed 
setting are made after this point. 

2. Turn rudder quickly by 15°, and keep it there until steady yaw rate is 
maintained for one minute. 

3. Reduce rudder angle by 5°, and keep it there until steady yaw rate is 
maintained for one minute. 

4. Repeat in decrements of -5°, to -15°. 

5. Proceed back up to 15°. 

The Dieudonne maneuver has the vessel path following a growing spiral, and 
then a contracting spiral in the opposite direction. The test reveals if the vessel 
has a memory effect, manifested as a hysteresis in yaw rate r. For example, 
suppose that the hrst 15° rudder deflection causes the vessel to turn right, but 
that the yaw rate at zero rudder, on the way negative, is still to the right. The 
vessel has gotten “stuck” here, and will require a negative rudder action to 
pull out of the turn. But if the corrective action causes the vessel to turn left 
at all, the same memory effect may occur. It is easy to see that the rudder in 
this case has to be used excessively driving the vessel back and forth. We say 
that the vessel is unstable, and clearly a poor design. 

7.2 Zig-Zag Maneuver 

1. With zero rudder, achieve steady speed for one minute. 

2. Deflect the rudder to 20°, and hold until the vessel turns 20°. 

3. Deflect the rudder to -20°, and hold until the vessel turns to -20o with 
respect to the starting heading. 



4. Repeat. 




7.3 Circle Maneuver 
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This maneuver establishes several important characteristics of the yaw re- 
sponse. These are: the response time (time to reach a given heading), the yaw 
overshoot (amount the vessel exceeds ±20° when the rudder has turned the 
other way), and the total period for the 20° oscillations. Of course, similar 
tests can be made with different rudder angles and different threshold vessel 
headings. 

7.3 Circle Maneuver 

From a steady speed, zero yaw rate condition, the rudder is moved to a new 
setting. The vessel responds by turning in a circle. After steady state is 
reached again, parameters of interest are the turning diameter, the drift angle 
/5, the speed loss, and the angle of heel ip. 

7.3.1 Drift Angle 

The drift angle is the equivalent to angle of attack for lifting surfaces, and 
describes how the vessel “skids” during a turn. If the turning circle has radius 
R (measured from the vessel origin), then the speed tangential to the circle 
is [/ = rR. The vessel-reference velocity components are thus u = U cos [3 
and V = —U sin/3. A line along the vessel centerline reaches closest to the 
true center of the turning circle at a point termed the turning center. At 
this location, which may or may not exist on the physical vessel, there is no 
apparent lateral velocity, and it appears to an observer there that the vessel 
turns about this point. 

7.3.2 Speed Loss 

Speed loss occurs primarily because of drag induced by the drift angle. A 
vessel which drifts very little may have very little speed loss. 

7.3.3 Heel Angle 

Heel during turning occurs as a result of the intrinsic coupling of sway, yaw, 
and roll caused by the center of gravity. In a surface vessel, the fluid forces act 
below the waterline, but the center of gravity is near the waterline or above. 
When the rudder is hrst deflected, inertial terms dominate (Phase 1) and the 
sway equation is 
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{m — Yi,)v — {Yr — raxcY' = Y^d. (89) 

The coefficients for r are quite small, and thus the vessel hrst rolls to starboard 
(positive) for a positive rudder action. 

When steady turning conditions are reached (Phase 3), hydrodynamic forces 
equalize the centrifugal force mU r and the rudder force Ys6. The sway equation 

is 



-Y^v + {mU - Yr)r = (90) 

with Yr small, u < 0 when r > 0 for most vessels, and \Yrv\ > |W(5|. Because 
the centrifugal force acts above the waterline, the vessel ultimately rolls to 
port (negative) for a positive rudder action. 

The transition between the inertially-dominated and steady-turning regimes 
includes an overshoot: in the above formulas, the vessel overshoots the hnal 
port roll angle as the vessel slows. From the sway equation, we see that if 
the rudder is straightened out at this point, the roll will momentarily be even 
worse! 

In summary, the vessel rolls into the turn initially, but then out of the turn in 
the steady state. 

7.3.4 Heeling in Submarines with Sails 

Submarines typically roll into a turn during all phases. Unlike surface vessels, 
which have the rigid mass center above the center of fluid forcing, submarines 
have the mass center below the rudder action point, and additionally feel the 
effects of a large sail above both. The inertial equation 



(m — Yi,)v — (Yr — rrixc)! = Ys6 (91) 

is dominated by mi) (acting low), —Yi,v (acting high), and Ys6 (intermediate). 
Because IWI > 'm, the vessel rolls under the sail, the keel out of the turn. In 
the steady state. 



-YrV + (mU - Yr)r = Ys5. (92) 

The drift angle f3 keeps the W-force, acting high, toward the center of the turn, 
and again centrifugal force mU r causes the bottom of the submarine to move 
out of the turn. Hence, the roll angle of a submarine with a sail is always into 
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the turn, both initially and in the steady state. The heel angle declines as the 
speed of the submarine drops. 
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8 STREAMLINED BODIES 

8.1 Nominal Drag Force 

A symmetric streamlined body at zero angle of attack experiences only a drag 
force, which has the form 



Fa = -\pCaAoU\ (93) 

The drag coefficient Ca has both pressure and skin friction components, and 
hence area A^ is usually that of the wetted surface. Note that the A-subscript 
will be used to denote zero angle of attack conditions; also, the sign of Fa is 
negative, because it opposes the vehicle’s x-axis. 

8.2 Munk Moment 

Any shape other than a sphere generates a moment when inclined in an invis- 
cid flow. d’Alembert’s paradox predicts zero net force, but not necessarily a 
zero moment. This Munk moment arises for a simple reason, the asymmetric 
location of the stagnation points, where pressure is highest on the front of 
the body (decelerating flow) and lowest on the back (accelerating flow). The 
Munk moment is always destabilizing, in the sense that it acts to turn the 
vehicle perpendicular to the flow. 

Consider a symmetric body with added mass components along the vehicle 
(slender) x-axis (forward), and A^z along the vehicle’s 2 ;-axis 2 : (up). We will 
limit the present discussion to the vertical plane, but similar arguments can 
be used to describe the horizontal plane. Let a represent the angle of attack, 
taken to be positive with the nose up - this equates to a negative pitch angle 
(j) in vehicle coordinates, if it is moving horizontally. The Munk moment is: 

Mm = -]^{Azz - Az^z^U"^ An2a (94) 

^ -{Azz- Az,x)U^a. 

Azz > Axx for a slender body, and the negative sign indicates a negative 
pitch with respect to the vehicle’s pitch axis. The added mass terms Azz and 
Axx can be estimated from analytical expressions (available only for regular 
shapes such as ellipsoids), from numerical calculation, or from slender body 
approximation (to follow). 




8.3 Separation Moment 
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8.3 Separation Moment 

In a viscous fluid, flow over a streamlined body is similar to that of potential 
flow, with the exceptions of the boundary layer, and a small region near the 
trailing end. In this latter area, a helical vortex may form and convect down- 
stream. Since vortices correlate with low pressure, the effect of such a vortex 
is stabilizing, but it also induces drag. The formation of the vortex depends on 
the angle of attack, and it may cover a larger area (increasing the stabilizing 
moment and drag) for a larger angle of attack. For a small angle of attack, the 
transverse force Fn can be written in the same form as for control surfaces: 

Fn = ]^pCnAoU‘^ (95) 

■ 2^ 9a ° • 

With a positive angle of attack, this force is in the positive ^-direction. The 
zero-a drag force Fa is modihed by the vortex shedding: 

Fa = -^pCaAoU^ , where (96) 

Ca = Ca cos2 0. 

The last relation is based on writing Ca{U cos (pY as (Ca cos^ 0)7/^, i.e., a 
decomposition using apparent velocity. 

8.4 Net Effects: Aerodynamic Center 

The Munk moment and the moment induced by separation are competing, 
and their magnitudes determine the stability of a hullform. First we simplify: 



Fa '^a 

Fn 

Tfrra ~ 



Each constant 7 is taken as positive, and the signs reflect orientation in the 
vehicle reference frame, with a nose-up angle of attack. The Munk moment is a 
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pure couple which does not depend on a reference point. We pick a temporary 
origin O for Fn however, and write the total pitch moment about O as: 



M — Mjn + Fnln (97) 

( dm ~h 

where denotes the (positive) distance between O and the application point 
of Fn- The net moment about O is zero if we select 

/n = — , (98) 

'In 

and the location of O is then called the aerodynamic center or AC. 

The point AC has an intuitive explanation: it is the location on the hull where 
Fn would act to create the total moment. Hence, if the vehicle’s origin lies in 
front of AC, the net moment is stabilizing. If the origin lies behind AC, the 
moment is destabilizing. For self-propelled vehicles, the mass center must be 
forward of AC for stability. Similarly, for towed vehicles, the towpoint must 
be located forward of AC . In many cases with very streamlined bodies, the 
aerodynamic center is signihcantly ahead of the nose, and in this case, a rigid 
sting would have to extend at least to AC in order for stable towing. As a 
hnal note, since the Munk moment persists even in inviscid flow, AC moves 
inhnitely far forward as viscosity effects diminish. 

8.5 Role of Fins in Moving the Aerodynamic Center 

Control surfaces or hxed fins are often attached to the stern of a slender vehicle 
to enhance directional stability. Fixed surfaces induce lift and drag on the 
body: 



L = ^pA,U^Q(a) 
D = -^pA;U^Ci 



~ '^LOi 

— —Id (constant) 



(99) 



These forces act somewhere on the £n, and are signed again to match the 
vehicle frame, with 7 > 0 and a > 0. The summed forces on the body are 
thus: 
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X = Fa — |-D| cos a + |L| sin q; (100) 

^ -7a - 7 d + 

Z = F„ + |F| COSO; + |F| sina 
- 7n« + iLOi + 7D«- 

All of the forces are pushing the vehicle up. If we say that the hns act a 
distance If behind the temporary origin O, and that the moment carried by 
the hns themselves is very small (compared to the moment induced by Llf) 
the total moment is as follows: 

M = {-yrn + lnln)oi + {y 1 + 10)1 f a. ( 101 ) 

The moment about O vanishes if 



Im Inin + IfilfL + 7o)' (102) 

The Munk moment im opposes the aggregate effects of vorticity lift in and 
the hns’ lift and drag 7 l + 7d- With very large hns, this latter term is large, so 
that If might be very small; this is the case of AC moving aft toward the hns. 
A vehicle with excessively large hns will be difficult to turn and maneuver. 
Equation 102 contains two length measurements, referenced to an arbitrary 
body point O. To solve it explicitly, let Ifn denote the (positive) distance that 
the hns are located behind F„; this is likely a small number, since both ehects 
usually act near the stern. We solve for lf\ 



lf = 



'y+a 'Jnlfn 



(103) 



In + lL + Id 

This is the distance that AC is located forward of the hns, and thus AC can 
be referenced to any other hxed point easily. Without hns, it can be recalled 
that 



r 'ym 

n • 

7n 



Hence, the hns act directly in the denominator to shorten If. Note that if 
the hns are located forward of the vortex shedding force F„, i.e., Ifn < 0, If is 
reduced, but since AC is referenced to the hns, there is no net gain in stability. 
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8.6 Aggregate Effects of Body and Fins 

Since all of the terms discussed so far have the same dependence on a, it is 
possible to group them into a condensed form. Setting Fa and to account 
for the fuselage and hns, we have 



X = Fa--^pCaAaU^ (104) 

^ = Fn-^pC'^AoU^a 

M = -FnXAc - -^C'^AoU^XAca- 

8.7 Coefficients My^, Zq, and Mq for a Slender Body 

The angle of attack a is related to the cross-body velocity w as follows: 



a = —tan ^ (105) 

~ — — tor t/ >> w. 

We can then write several linear hydrodynamic coefficients easily: 

Zw = -^pC'aAoU (106) 

My, = ^pC'^AoUxAc- 

The rotation of the vessel involves complex flow, which depends on both w and 
g, as well as their derivatives. To start, we consider the contribution of the hns 
only - slender body theory, introduced shortly, provides good results for the 
hull. The hn center of pressure is located a distance If aft of the body origin, 
and we assume that the vehicle is moving horizontally, with an instantaneous 
pitch angle of 6. The angle of attack seen by the hn is a combination of a part 
due to 9 and a part linear with q: 



a,~-e+’XL 



( 107 ) 




8.7 CoefEcients Zy^, My,, Zq, and Mq for a Slender Body 
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and so lateral force and moment derivatives (for the fin alone) emerge as 

z, = -^-pc;a,ui, ( 108 ) 

M, = -\pC[A,UP,. 
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9 SLENDER-BODY THEORY 

9.1 Introduction 

Consider a slender body with d « L, that is mostly straight. The body 
could be asymmetric in cross-section, or even flexible, but we require that 
the lateral variations are small and smooth along the length. The idea of 
the slender-body theory, under these assumptions, is to think of the body as 
a longitudinal stack of thin sections, each having an easily-computed added 
mass. The effects are integrated along the length to approximate lift force and 
moment. Slender-body theory is accurate for small ratios d/L, except near 
the ends of the body. 

As one example, if the diameter of a body of revolution is d{s), then we can 
compute Sma{x), where the nominal added mass value for a cylinder is 

dnia = p—d^5x. (109) 

The added mass is equal to the mass of the water displaced by the cylinder. 
The equation above turns out to be a good approximation for a number of 
two-dimensional shapes, including flat plates and ellipses, if d is taken as the 
width dimension presented to the flow. Many formulas for added mass of two- 
dimensional sections, as well as for simple three-dimensional bodies, can be 
found in the books by Newman and Blevins. 

9.2 Kinematics Following the Fluid 

The added mass forces and moments derive from accelerations that a fluid 
particles experience when they encounter the body. We use the notion of a 
fluid derivative for this purpose: the operator d/dt indicates a derivative taken 
in the frame of the passing particle, not the vehicle. Hence, this usage has an 
indirect connection with the derivative described in our previous discussion of 
rigid-body dynamics. 

For the purposes of explaining the theory, we will consider the two-dimensional 
heave/surge problem only. The local geometry is described by the location of 
the centerline; it has vertical location (in body coordinates) of Zb{x,t), and 
local angle a{x,t). The time-dependence indicates that the configuration is 
free to change with time, i.e., the body is flexible. Note that the curvilinear 
coordinate s is nearly equal to the body-reference (linear) coordinate x. 

The velocity of a fluid particle normal to the body at x is Wn(t, x): 




9.3 Derivative Following the Fluid 
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w 



n 



dzb . 

— — cos a — U sm a. 
dt 



( 110 ) 



The first component is the time derivative in the body frame, and the second 
due to the deflection of the particle by the inclined body. If the body reference 
frame is rotated to the flow, that is, if tc 7 ^ 0, then dzb/dt will contain w. For 
small angles, sin a ~ tana = dzh/dx, and we can write 



w 



n 






dzb _ ^jdzb 

dt dx 



The fluid derivative operator in action is as follows: 



w 



n 



dzb 

dt 





Zb- 



9.3 Derivative Following the Fluid 

A more formal derivation for the fluid derivative operator is quite simple. Let 
jj,{x, t) represent some property of a fluid particle. 



lim — [/i(a; + 5x, t + 5t) — fi(x, t)] 
dt 

^ 

dt dx 

The second equality can be verihed using a Taylor series expansion of fr{x + 
Sx, t + 6t)\ 

/ c r \ / \ c c 7 

fiix + ox, t + ot) = ^[x, t) + -^ot + —ox + h.o.t., 

dt dx 

and noting that Sx = —Ft St. The fluid is convected downstream with respect 
to the body. 
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9.4 Differential Force on the Body 

If the local transverse velocity is Wn(x, t), then the differential inertial force on 
the body here is the derivative (following the flnid) of the momentum; 

SF = -^[ma{x,t)wn{x,t)]5x. (Ill) 

Note that we could here let the added mass vary with time also - this is the 
case of a changing cross-section! The lateral velocity of the point zi,{x) in the 
body-reference frame is 



such that 



dzb 

dt 



= w — xq, 



Wn = w — xq — Ua. 
Taking the derivative, we have 



( 112 ) 



( 113 ) 



T = _ ("4 - U^ |m„(i, ()(«)(«) -xq(t) - Ua(x,t))] . 

We now restrict ourselves to a rigid body, so that neither rua nor a may change 
with time. 



Sx 



ma{x){—w + xq) -|- U 



A 

dx 



ma{x){w 



xq — Ua)] . 



(114) 



9.5 Total Force on a Vessel 

The net lift force on the body, computed with strip theory is 



fXN 

z= SFdx (115) 

J Xjp 

where xt represents the coordinate of the tail, and a ; at is the coordinate of the 
nose. Expanding, we have 




9.6 Total Moment on a Vessel 
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rXN rXN Q 

Z= ma{x) [-W + xq] dx + U / — [ma{x){w — xq — Ua)] dx 

Jx'T J xjp UX 

= -m33W - ms^q + Uma{x){w - xq - Ua)\lZl^. 

We made use here of the added mass dehnitions 



rXN 

mss = / ma{x)dx 

J X'T 

PXn 

77135 = — xma{x)dx. 

J XT 

Additionally, for vessels with pointed noses and flat tails, the added mass ma 
at the nose is zero, so that a simpler form occurs: 

Z = —mssw — ms5q — Uma{xT){w — xxq — Ua{xT)). (116) 
In terms of the linear hydrodynamic derivatives, the strip theory thus provides 



Zw ^Xlss 



Zq — —77735 

Zw = -Umaixr) 

Zq = UxrmaixT) 
Za(xT) = U‘^ma{xT). 



It is interesting to note that both Zw and Za(xT) depend on a nonzero base 
area. In general, however, potential flow estimates do not create lift (or drag) 
forces for a smooth body, so this should come as no surprise. The two terms are 
clearly related, since their difference depends only on how the body coordinate 
system is oriented to the flow. Another noteworthy fact is that the lift force 
depends only on a at the tail; a could take any value (s) along the body, with 
no effect on Z . 



9.6 Total Moment on a Vessel 

A similar procedure can be applied to the moment predictions from slender 
body theory (again for small a): 
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M = - 



x 6 Fdx 



r^N / (9 d \ 

I X \ — - U—j [ma{x){w - xq - Ua)]dx 



rXN rXN Q 

/ xma{x){w — xq)dx — U / x— [ma{x){w — xq + Ua)] dx 

Jx'T J UX 



Then we make the further dehnition 



= 



x‘^ma{x)dx, 



(note that mss = ^53) and use integration by parts to obtain 

M = -m^^w - m^^q - U xnia{x){w - xq - U + 

rXN 

U / ma{x){w — xq — Ua)dx. 

J XT' 

The integral above contains the product ma{x)a{x), which must be calculated 
if a changes along the length. For simplicity, we now assume that a is in fact 
constant on the length, leading to 



M = —mz^w — m^^q + UxTma{xT){w — XTq — Ua) + 

Um^sw + Um^^q - U'^mssa. 

Finally, the linear hydrodynamic moment derivatives are 



= -mss 
Mq = -mss 

Mu, = UxTfriaixT) + Umss 
Mq = —Ux'^ma{xT) + Urri35 
Ma = -U‘^XTma{xT) ~ 

The derivative M^, is closely-related to the Munk moment, whose linearization 
would provide M^, = (mss — rnu)U. The Munk moment (an exact result) may 
therefore be used to make a correction to the second term in the slender-body 
approximation above of M^,. As with the lift force, M^, and are closely 
related, depending only on the orientation of the body frame to the flow. 




9. 7 Relation to Wing Lift 
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9.7 Relation to Wing Lift 

There is an important connection between the slender body theory terms in- 
volving added mass at the tail {ma{xT)), and low aspect-ratio wing theory. 
The lift force from the latter is of the form L = —^pLfAC[w, where A = cs, 
the product of chord (long) and span (short). 

The lift coefficient slope is approximated by (Hoerner) 

c; « tir(Afi), (117) 

where {AR) is the aspect ratio. Inserting this approximation into the lift 
formula, we obtain 

L = —^ps'^Uw. ( 118 ) 

Now we look at a slender body approximation of the same force: The added 
mass at the tail is nia{xT) = ps'^n/A, and using the slender-body estimate for 
Zu,, we calculate for lift: 



Z 



—ma{xT)Uw 

— ^ps^Uw. 

4 



Slender-body theory is thus able to recover exactly the lift of a low-aspect 
ratio wing. Where does the slender-body predict the force will act? Recalling 
that Mu, = 77m33 -|- and since 77133 = 0 for a front-back symmetric 

wing, the estimated lift force acts at the trailing edge. This location will tend 
to stabilize the wing, in the sense that it acts to orient the wing parallel to 
the incoming flow. 



9.8 Convention: Hydrodynamic Mass Matrix A 

Hydrodynamic derivatives that depend on accelerations are often written as 
components of a mass matrix A. By listing the body-referenced velocities 
in the order s = [u, v, w, p, g, r], we write (M -|- H)T = F, where M is the 
mass matrix of the material vessel and F is a generalized force. Therefore 
H33 = -Z^, ^5^3 = -M^, and so on. 
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10 PRACTICAL LIFT CALCULATIONS 

10.1 Characteristics of Lift-Producing Mechanisms 

At a small angle of attack, a slender body experiences transverse force due to: 
helical body vortices, the blunt trailing end, and hns. The helical body vortices 
are stable and symmetric in this condition, and are convected continuously into 
the wake. The low pressure associated with the vortices provides the suction 
force, usually toward the stern of the vehicle. The blunt trailing end induces 
lift as a product of added mass effects, and can be accurately modeled with 
slender body theory. A blunt trailing edge also induces some drag, which itself 
is stabilizing. The hns can often be properly modeled using experimental data 
parametrized with aspect ratio and several other geometric quantities. 

As the angle of attack becomes larger, the approximations in the hn and 
slender-body analysis will break down. The helical vortices can become bigger 
while remaining stable, but eventually will split randomly. Some of it convects 
downstream, and the rest peels away from the body; this shedding is nonsym- 
metric, and greatly increases drag by widening the wake. In the limit of a 90° 
angle of attack, vorticity sheds as if from a bluff-body, and there is little axial 
convection. 



10.2 Jorgensen’s Formulas 

There are some heuristic and theoretical formulas for predicting transverse 
force and moment on a body at various angles of attack, and we now present 
one of them due to Jorgensen. The formulas provide a good systematic proce- 
dure for design, and are best suited to vehicles with a blunt trailing edge. We 
call the area of the stern the base area. 

Let the body have length L, and reference area Aj.. This area could be the 
frontal projected area, the planform area, or the wetted area. The body travels 
at speed U, and angle of attack a. The normal force, axial force, and moment 
coefficients are dehned as follows: 



Cn 

Ca 



Fn 

\pU^Ar 

Fa 

IpU^Ar 



( 119 ) 
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‘^m 

^ ~ IpU^ArL' 

The moment is taken about a point Xm, measured back from the nose; this 
location is arbitrary, and appears explicitly in the formula for Cm- Jorgensen 
gives the coefficients as follows: 



Cn 

Ca 

Cm 



Ab . „ O' Ap . 2 

— smiacos — \- -ACd sm a 
4 9 4 

Cao cos^ a 

V — Ab{L — Xr, 



ArL 



sm 2a cos 2 ~ a 






(120) 
( 121 ) 
sin^ q;.(122) 



We have listed only the formulas for the special case of circular cross-section, 
although the complete equations do account for more complex shapes. Further, 
we have assumed that L » D, which is also not a constraint in the complete 
equations. The parameters used here are 



• Ab- stern base area. Ab = 0 for a body that tapers to a point at the 
stern. 

• Cd„'- crossflow drag coefficient; equivalent to that of an inhnite circular 
cylinder. If “normal” Reynolds number i?e„ = U sin aD/u, 

- Cd- 1.2, Ren < 3 X 10^ 

- Cd^ 0.3, 3 X 10^ < Rcn < 7 X 10® 

- Cd^ 0.6, 7 X 10^ < Ren- 

• Ap\ planform area. 

• Cao'- axial drag at zero angle of attack, both frictional and form. Cao — 
0.002 — 0.006 for slender streamlined bodies, based on wetted surface 
area. It depends on Re = UL/u. 

• V; body volume. 

• Xc'- distance from the nose backwards to the center of the planform area. 
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In the formula for normal force, we see that if Ai, = 0, only drag forces act 
to create lift, through a sin^ a-dependence. Similarly, the axial force is simply 
the zero-a result, modihed by cos^a. In both cases, scaling of into the 
body principle directions is all that is required. 

There are several terms that match exactly the slender-body theory approxi- 
mations for small a. These are the hrst term in the normal force (Cat), and 
the entire hrst term in the moment {Cm), whether or not Ab = 0. Finally, we 
note that the second term in Cm disappears if Xm = Xc, i.e., if the moment is 
referenced to the center of the planform area. The idea here is that the fore 
and aft components of crosshow drag cancel out. 

The aerodynamic center (again referenced toward the stern, from the nose) 
can be found after the coefficients are computed: 

xac = Xm + (123) 

Oat 

As written, the moment coefficient is negative if the moment destabilizes the 
body, while C^ is always positive. Thus, the moment seeks to move the AC 
forward on the body, but the effect is moderated by the lift force. 

10.3 Hoerner’s Data: Notation 

An excellent reference for experimental data is the two- volume set by S. Ho- 
erner. It contains a large amount of aerodynamic data from many different 
types of vehicles, wings, and other common engineering shapes. A few nota- 
tions are used throughout the books, and are described here. 

First, dynamic pressure is given as g = such that two typical body lift 

coefficients are: 



F 

I)Lq 
Y 

Wq' 

The first version uses the rectangular planform area as a reference, while the 
second uses the square frontal area. Hence, Cy = Cy^D/L. Two moment 
coefficients are: 



Cy = 




10.4 Slender-Body Theory vs. Experiment 
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r - ^ 

c 

LD^q' 

where N is the moment taken about the body mid-point, and Ni is taken 
about the nose. Note that the reference area for moment is the square frontal 
projection, and the reference length is body length L. The following relation 
holds for these dehnitions 



Cni = Cn + -Cva- 

The lift and moment coefficients are strongly dependent on angle of attack; 
Hoerner uses the notation 



dCjy 
db 

OCm, 
db 

dCy 

db ’ 

and so on, where b is the angle of attack, usually in degrees. It follows from 
above that C^-b = Cnb + Cydb/‘2^- 

10.4 Slender-Body Theory vs. Experiment 

In an experiment, the net moment is measured, comprising both the desta- 
bilizing part due to the potential flow, and the stabilizing part due to vortex 
shedding or a blunt tail. Comparison of the measurements and the theory 
allows us to place the action point of the suction force. This section gives the 
formula for this location in Hoerner’s notation, and gives two further examples 
of how well the slender-body theory matches experiments. 

For L/D > 6, the slender-body (pure added mass) estimates give Cnb — 
— O.OlS/degf, acting to destabilize the vehicle. The value compares well with 
—0.027 /deg for a long cylinder and —O.OIS/ deg for a long ellipsoid; it also 



Cnb — 
Cn-b = 

Cyb = 




72 



10 PRACTICAL LIFT CALCULATIONS 



reduces to —0.009 for L/ D = 4. Note that the negative sign here is consistent 
with Hoerner’s convention that destabilizing moments have negative sign. 

The experimental lift force is typically given by Cyb — O.OOS/degf; this acts at 
a point on the latter half of the vehicle, stabilizing the angle. Because this 
coefficient scales roughly with wetted area, proportional to LD, it changes 
little with LID. It can be compared with a low-aspect ratio wing, which 
achieves an equivalent lift of ti{AR)/2 = 0.0027 for {AR) = 10 ~ D/L. 

The point at which the viscous forces act can then be estimated as the following 
distance aft of the nose: 



X 

I 



Cn-h — C, 



nb 



a 



ydb 



( 124 ) 



The calculation uses experimental values of Cydi, and the moment slope 
referenced to the nose. In the table following are values from Hoerner (p. 13-2, 
Figure 2) for a symmetric and a blunt-ended body. 





symmetric 


blunt 


L/D 


6.7 


6.7 


Cnb 


-0.012 


-0.012 


Cyb 


0.0031 


0.0037 


Cydb 


0.021 


0.025 


Cn-b 


0.0012 (stable) 


0.0031 (stable) 


Cnb 


-0.0093 (unstable) 


-0.0094 (unstable) 


xjL 


0.63 


0.60 



In comparing the two body shapes, we see that the moment at the nose is 
much more stable (positive) for the body with a blunt trailing edge. At the 
body midpoint, however, both vehicles are equally unstable. The blunt-tailed 
geometry has a much larger lift force, but it acts too close to the midpoint to 
add any stability there. 

The lift force dependence on the blunt tail is not difficult to see, using slender- 
body theory. Consider a body, with trailing edge radius r. The slender-body 
lift force associated with this end is simply the product of speed U and local 
added mass ma{xT) (in our previous notation). It comes out to be 



Z=^pU\7rr^){2a), (125) 

such that the hrst term in parentheses is an effective area, and the second is 
a lift coefficient. With respect to the area Trr^, the lift curve slope is therefore 
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2/ rad. Expressed in terms of the Hoerner reference area the equivalent 
lift coefficient is Cydb = O.OOdd/degf, where we made the assumption here that 
2r/D = 0.4 for the data in the table. This lift difference, due solely to the 
blunt end condition, is consistent with measurements. 

10.5 Slender-Body Approximation for Fin Lift 

Let us now consider two hns of span s each, acting at the tail end of the 
vehicle. This is the case if the vehicle body tapers to a point where the fins 
have their trailing edge. The slender-body approximation of lift as a result of 
blunt-end conditions is 



Z = xs^pU'^a. (126) 

Letting the aspect ratio be {AR) = [2s)‘^/Af, where Af is the total area of the 
£n pair, substitution will give a lift curve slope of 

C'i = l{AR). 

This is known as Jones’ formula, and is quite accurate for {AR) ~ 1. It is 
inadequate for higher-aspect ratio wings however, overestimating the lift by 
about 30% when {AR) = 2, and worsening further as {AR) grows. 
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11 FINS AND LIFTING SURFACES 

Vessels traveling at significant speed typically use rudders, elevators, and other 
streamlined control surfaces to maneuver. Their utility arises mainly from the 
high lift forces they can develop, with little drag penalty. Lift is always dehned 
to act in a direction perpendicular to the flow, and drag in the same direction 
as the flow. 

11.1 Origin of Lift 

A lifting surface is nominally an extrusion of a streamlined cross-section: the 
cross-section has a rounded leading edge, sharp trailing edge, and a smooth 
surface. The theory of lifting surfaces centers on the Kutta condition, which 
requires that fluid particle streamlines do not wrap around the trailing edge 
of the surface, but instead rejoin with streamlines from the other side of the 
wing at the trailing edge. This fact is true for a non-stalled surface at any 
angle of attack. 

Since the separation point on the front of the section rotates with the angle 
of attack, it is clear that the fluid must travel faster over one side of the 
surface than the other. The reduced Bernoulli pressure this induces can then 
be thought of as the lift-producing mechanism. More formally, lift arises from 
circulation T: 



T = j^V-ds. (127) 

and then L = —pUT. Circulation is the integral of velocity around the cross- 
section, and a lifting surface requires circulation in order to meet the Kutta 
condition. 

11.2 Three-Dimensional Effects: Finite Length 

Since all practical lifting surfaces have hnite length, the flow near the ends 
may be three-dimensional. Prandtl’s inviscid theory provides some insight. 
Since bound circulation cannot end abruptly at the wing end, it continues on 
in the fluid, leading to so-called wing-tip vortices. This continuation causes 
induced velocities at the tips, and some induced drag. Another description for 
the wing-tip vortices is that the pressure difference across the surface simply 
causes flow around the end. 
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A critical parameter which governs the extent of three-dimensional effects is 
the aspect ratio: 



span span‘s 
chord area 

The second representation is useful for non-rectangular control surfaces. The 
effective span is taken to be the length between the free ends of a symmetric 
wing. If the wing is attached to a wall, the effective span is twice the physical 
value, by reflection, and in this case the effective aspect ratio is therefore twice 
the physical value. 

The aspect ratio is a strong determinant of wing performance: for a given 
angle of attack, a larger aspect ratio achieves a higher lift value, but also stalls 
earlier. 

Lift is written as 




L = ^pU^ACi, (129) 

where A is the single-side area of the surface. For angles of attack a below 
stall, the lift coefficient Ci is nearly linear with a: Ci = C[a, where C'l is called 
the lift coefficient slope, and has one empirical description 

= 1 . d . 1 — ■ (130) 

27TO tt{AR) ' 2tt{ARP 

where a is in radians, a ~ 0.90, and AR is the effective aspect ratio. When 
AR — oo, the theoretical and maximum value for C'l is 27 t. 

The lift generated on a surface is the result of a distributed pressure held; this 
fact creates both a net force and a net moment. A single equivalent force acts 
at a so-called center of action xa, which depends mainly on the aspect ratio. 
For high AR, xa — c/4, measured back from the front of the wing. For low 
AR, Xa — cl2. 

11.3 Ring Fins 

Ring hns are useful when space allows, since they are omnidirectional, and 
structurally more robust than cantilevered plane surfaces. The effective aspect 
ratio for a ring of diameter d is given as 

arAA, 

T^C 



( 131 ) 
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The effective area of the ring is taken as 

. ^ 7 

-Ag —dc^ 

and we thus have L = ^pU'^AgCla, where one formula for C'l is 



C',= 



0.63 + 



r{AR) 



(132) 



(133) 
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12 PROPELLERS AND PROPULSION 

12.1 Introduction 

We discuss in this section the nature of steady and unsteady propulsion. In 
many marine vessels and vehicles, an engine (diesel or gas turbine, say) or an 
electric motor drives the propeller through a linkage of shafts, reducers, and 
bearings, and the effects of each part are important in the response of the net 
system. Large, commercial surface vessels spend the vast majority of their time 
operating in open-water and at constant speed. In this case, steady propulsion 
conditions are generally optimized for fuel efficiency. An approximation of the 
transient behavior of a system can be made using the quasi-static assumption. 
In the second section, we list several low-order models of thrusters, which have 
recently been used to model and simulate truly unsteady conditions. 

12.2 Steady Propulsion of Vessels 

The notation we will use is as given in Table 1, and there are two different flow 
conditions to consider. Self-propelled conditions refer to the propeller being 
installed and its propelling the vessel; there are no additional forces or moments 
on the vessel, such as would be caused by a towing bar or hawser. Furthermore, 
the flow around the hull interacts with the flow through the propeller. We use 
an sp subscript to indicate specihcally self-propulsion conditions. Conversely, 
when the propeller is run in open water, i.e., not behind a hull, we use an o 
subscript; when the hull is towed with no propeller we use a t subscript. When 
subscripts are not used, generalization to either condition is implied. Finally, 
because of similitude (using diameter D in place of L when the propeller is 
involved), we do not distinguish between the magnitude of forces in model and 
full-scale vessels. 

12.2.1 Basic Characteristics 

In the steady state, force balance in self-propulsion requires that 

Rsp — Tgp. (134) 

The gear ratio A is usually large, indicating that the propeller turns much 
more slowly than the driving engine or motor. The following relations dehne 
the gearbox: 
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Tgp 


N 


hull resistance under self-propulsion 


Rt 


N 


towed hull resistance (no propeller attached) 


T 


N 


thrust of the propeller 


Ue 


Hz 


rotational speed of the engine 


'l^m 


Hz 


maximum value of Ug 


Up 


Hz 


rotational speed of the propeller 


X 




gear ratio 


Qe 


Nm 


engine torque 


Qp 


Nm 


propeller torque 


Vg 




gearbox efficiency 


Pe 


W 


engine power 


Pp 


w 


propeller shaft power 


D 


m 


propeller diameter 


U 


m/s 


vessel speed 


Up 


m/s 


water speed seen at the propeller 


Qm 


Nm 


maximum engine torque 


f 


kg/s 


fuel rate (or energy rate in electric motor) 


fm 


kg/s 


maximum value of / 



Table 1: Nomenclature 



He = \Up (135) 

Qp 

and power follows as Pp = rjgPg, for any flow condition. We call J = Up/upD 
the advance ratio of the prop when it is exposed to a water speed Up] note 
that in the wake of the vessel, Up may not be the same as the speed of the 
vessel U. A propeller operating in open water can be characterized by two 
nondimensional parameters which are both functions of J: 



Kt = 



Kq = 



To 

Qpo 

pnlD^ 



(thrust coefficient) 



(136) 



(torque coefficient). 



(137) 
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Figure 4: Typical thrust and torque coefficients. 



The open propeller efficiency can be written then as 



TM _ .J{U)Kt 

27TTlpC^p^ 2 ttKq 



(138) 



This efficiency divides the useful thrust power by the shaft power. Thrust and 
torque coefficients are typically nearly linear over a range of J, and therefore 
£t the approximate form: 



Kt{J) = P 1 -P 2 J (139) 

^q(J) = 71 - 72 ^- 

As written, the four coefficients [Pi, P 2 , 7i, 72] are usually positive, as shown in 
the figure. 

We next introduce three factors useful for scaling and parameterizing our 
mathematical models: 

• Up = U{1 — w)] tc is referred to as the wake fraction. A typical wake 
fraction of 0.1, for example, indicates that the incoming velocity seen by 
the propeller is only 90% of the vessel’s speed. The propeller is operating 
in a wake. 

In practical terms, the wake fraction comes about this way: Suppose the 
open water thrust of a propeller is known at a given U and Up. Behind 
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a vessel moving at speed U, and with the propeller spinning at the same 
Up, the prop creates some extra thrust, w scales U at the prop and thus 
J] w is then chosen so that the open water thrust coefficient matches 
what is observed. The wake fraction can also be estimated by making 
direct velocity measurements behind the hull, with no propeller. 

• Rt = Rspil — t). Often, a propeller will increase the resistance of the 
vessel by creating low-pressure on its intake side (near the hull), which 
makes Rsp > Rt- In this case, t is a small positive number, with 0.2 as a 
typical value, t is called the thrust deduction even though it is used to 
model resistance of the hull; it is obviously specihc to both the hull and 
the propeller(s), and how they interact. 

The thrust deduction is particularly useful, and can be estimated from 
published values, if only the towed resistance of a hull is known. 

• Qpo = VrQpsp- The rotative efhciency rjR, which may be greater than 
one, translates self-propelled torque to open water torque, for the same 
incident velocity Up, thrust T, and rotation rate rtp. rjn is meant to 
account for spatial variations in the wake of the vessel which are not 
captured by the wake fraction, as well as the turbulence induced by the 
hull. Note that in comparison with the wake fraction, rotative efficiency 
equalizes torque instead of thrust. 

A common measure of efficiency, the quasi-propulsive efficiency, is based on 
the towed resistance, and the self-propelled torque. 



Vqp 



RtU 



2TTnpQp^^ 

Toil- t)UpT]R 



27mp(l - w)Qp^ 

(l-t) 



VoVr 



(l-w)' 



(140) 



To and Qp^ are values for the inflow speed Up, and thus that rjo is the open 
propeller efficiency at this speed. It follows that To{Up) = Tgp, which was used 
to complete the above equation. The quasi-propulsive efficiency can be greater 
than one, since it relies on the towed resistance and in general Rt > Rsp- The 
ratio (1 — t)/(l — w) is often called the hull efficiency, and we see that a 
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small thrust deduction t and a large wake fraction w are beneficial effects, but 
which are in competition. A high rotative efficiency and open water propeller 
efficiency (at Up) obviously contribute to an efficient overall system. 

12.2.2 Solution for Steady Conditions 

The linear form of and Kq (Equation 140) allows a closed-form solntion 
for the steady-operating conditions. Suppose that the towed resistance is of 
the form 



Rt — -pCrA^U"^, (141) 

where Cr is the resistance coefficient (which will generally depend on Re and 
Fr), and is the wetted area. Equating the self-propelled thrust and resis- 
tance then gives 



T 

sp 

T 



KT{ J(UypnlD‘{l - 1) 



(ft - ftJ(C/p))pn|C‘(l - t) 



ft - ft J(ftp) 



J{U,) 



Rsp 

l^pCrA^U"^ 



f/2 

P 



pCrA^u^ 

2 (1 — w)^ 



C A 

^r-^w 



2F)2(1 -t)(l -u;)2 



J(Up)^ 



-^2 + \Z/d| + 4/?l5 



26 



(142) 



The last equation predicts the steady-state advance ratio of the vessel, depend- 
ing only on the propeller open characteristics, and on the hnll. The vessel speed 
can be computed by recalling that J{U) = U /upD and Up = U {1 — w), bnt it 
is clear that we need now to hnd Up. This requires a torque equation, which 
necessitates a model of the drive engine or motor. 



12.2.3 Engine/Motor Models 

The torque-speed maps of many engines and motors £t the form 




82 



12 PROPELLERS AND PROPULSION 




Figure 5: Typical gas turbine engine torque-speed characteristic for increasing 
fuel rates fi,f 2 ,fs- 



Qe = QmF{f/fm, rie/n^), (143) 

where F() is the characteristic function. For example, gas turbines roughly £t 
the curves shown in the hgure (Rubis). More specifically, if Fi) has the form 



F{f/ffm,ne/nm) 




He 

— + 
'^m 



— ttl 



Up 



a2- 




( 144 ) 



then a closed-form solution for Ue (and thus Up) can be found. The manipula- 
tions begin by equating the engine and propeller torque: 



QpoiJiUp)) 

pnlD^KgiJiUp)) 

pnlD^{^,-^2J{Up)) 

nl 

(’^m/A)2 



A) 



VRQpspi'^i^p)) 

VRVg^Qe 

VRVg^QmF{f/ fm, ne/rim) 

VRVg^Qm f Rp 

pB5(7i - l^J(U,))(nJXP 

N/* 

€ 

— €CHi p 6^0;^ -|- 46C^ 



2 



( 145 ) 
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Note that the fuel rate enters through both ai and 0 : 2 . 

The dynamic response of the coupled propulsion and ship systems, under the 
assumption of quasi-static propeller conditions, is given by 



(m + ma)u = Tsp - Rsp (146) 

^Tllpflp Pgl^Qe QpsP' 

Making the necessary substitutions creates a nonlinear model with / as the 
input; this is left as a problem for the reader. 

12.3 Unsteady Propulsion Models 

When accurate positioning of the vehicle is critical, the quasi-static assumption 
used above does not suffice. Instead, the transient behavior of the propulsion 
system needs to be considered. The problem of unsteady propulsion is still in 
development, although there have been some very successful models in recent 
years. It should be pointed out that the models described below all pertain to 
open-water conditions and electric motors, since the positioning problem has 
been central to bluff vehicles with multiple electric thrusters. 

We use the subscript m to denote a quantity in the motor, and p for the 
propeller. 

12.3.1 One-State Model: Yoerger et ai. 

The torque equation at the propeller and the thrust relation are 



IpUJp 

T = CtUp\up\. (148) 

where Ip is the total (material plus fluid) inertia reflected to the prop;the 
propeller spins at ujp radians per second. The differential equation in Up pits 
the torque delivered by the motor against a quadratic-drag type loss which 
depends on rotation speed. The thrust is then given as a static map directly 
from the rotation speed. 

This model requires the identihcation of three parameters: Ip, K,_^, and Ct- 
It is a hrst-order, nonlinear, low-pass hlter from Qm to T, whose bandwidth 
depends directly on Qm- 




84 



12 PROPELLERS AND PROPULSION 



12.3.2 Two-State Model: Healey et al. 

The two-state model includes the velocity of a mass of water moving in the 
vicinity of the blades. It can accommodate a tunnel around the propeller, 
which is very common in thrusters for positioning. The torque equation, sim- 
ilarly to the above, is referenced to the motor and given as 

^mOrn T Qp/ (149) 

Here, represents losses in the motor due to spinning (friction and resistive), 
and Ky is the gain on the input voltage (so that the current ampliher is included 
in Ky). The second dynamic equation is for the fluid velocity at the propeller: 

pAL-filp = -pAA(3{Up -U)\Up-U\+T. (150) 

Here A is the disc area of the tunnel, or the propeller disk diameter if no tunnel 
exists. L denotes the length of the tunnel, and 7 is the effective added mass 
ratio. Together, pALj is the added mass that is accelerated by the blades; this 
mass is always nonzero, even if there is no tunnel. The parameter A /3 is called 
the differential momentum flux coefficient across the propeller; it may be on 
the order of 0.2 for propellers with tunnels, and up to 2.0 for open propellers. 
The thrust and torque of the propeller are approximated using wing theory, 
which invokes lift and drag coefficients, as well as an effective angle of attack 
and the propeller pitch. However, these formulae are static maps, and therefore 
introduce no new dynamics. As with the one-state model of Yoerger et al, this 
version requires the identihcation of the various coefficients from experiments. 
This model has the advantage that it creates a thrust overshoot for a step 
input, which is in fact observed in experiments. 
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Vehicles which are towed have some similarities to the vehicles that have been 
discussed so far. For example, towed vehicles are often streamlined, and usu- 
ally need good directional stability. Some towed vehicles might have active 
lifting surfaces or thrusters for attitude control. On the other hand, if they 
are to be supported by a cable, towed vehicles may be quite heavy in water, 
and do not have to be self-propelled. The cable itself is an important factor in 
the behavior of the complete towed system, and in this section, we concentrate 
on cable mechanics more than vehicle characteristics, which can generally be 
handled with the same tools as other vehicles, i.e., slender-body theory, wing 
theory, linearization, etc.. Some basic guidelines for vehicle design are given 
at the end of this section. 

Modern cables can easily exceed 5000m in length, even a heavy steel cable with 
2cm diameter. The cables are generally circular in cross section, and may 
carry power conductors and multiple communication channels (fiber optic). 
The extreme L/ D ratio for these cables obviates any bending stiffness effects. 
Cable systems come in a variety of configurations, and one main division may 
be made simply of the density of the cable. Light-tether systems are character- 
ized by neutrally-buoyant (or nearly so) cables, with either a minimal vehicle 
at the end, as in a towed array, or a vehicle capable of maneuvering itself, such 
as a remotely-operated vehicle. The towed array is a relatively high-velocity 
system that nominally streams out horizontally behind the vessel. An ROV, 
on the other hand, operates at low speed, and must have large propulsors to 
control the tether if there are currents. Heavy systems, in contrast, employ a 
heavy cable and possibly a heavy weight; the rationale is that gravity will tend 
to keep the cable vertical and make the deployment robust against currents 
and towing speed. The heavy systems will generally transmit surface motions 
and tensions to the towed vehicle much more easily than light-tether systems. 
We will not discuss light systems specifically here, but rather look at heavy 
systems. Most of the analysis can be adapted to either case, however. 

13.1 Statics 

13.1.1 Force Balance 

For the purposes of deriving the static configuration of a cable in a flow, we 
assume for the moment that that it is inextensible. Tension and hydrostatic 
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pressure will elongate a cable, but the effect is usually a small percentage of 
the total length. 

We employ the curvilinear axial coordinate s, which we take to be zero at the 
bottom end of the cable; upwards along the cable is the positive direction. 
The free-body diagram shown has the following components: 

• Wn'- net in- water weight of the cable per unit length. 

• Rn{s): external normal force, per unit length. 

• Rt{s): external tangential force, per unit length. 

• T{s): local tension. 

• 4>{s): local inclination angle. 

Force balance in the tangential and normal coordinates gives two coupled 



amp; = amp; Wn sin (j) — Rt 

amp] = amp] Wn cos 0 + 

The external forces are primarily fluid drag; the tangential drag is controlled 
by a frictional drag coefficient Ct, and the normal drag scales with a crossflow 
drag coefficient In both cases, the fluid velocity vector, U horizontal 
toward the left, is to be projected onto the relevant axes, leading to 



(151) 

(152) 



equations for T and <p: 

dT 

ds 

ds 
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■ —^pCtdU'^ cos^ 0 


(153) 


— ^pCndU'^ sin^ 0. 


(154) 



Note that we simplified the drag laws from a usual form u|u| to since as 
drawn, 0 < 0 < tt/2. 

The equations for T and 0 can be integrated along the cable coordinate s to 
find the cable’s static configuration. Two boundary conditions are needed, and 
the common case is that a force balance on the vehicle, dominated by drag, 
weight, and the cable tension, provides both T(0) and 0(0). For example, a 
very heavy but low-drag vehicle will impose 0(0) ~ tt/ 2, with T(0) equal to 
the in-water weight of the vehicle. 

With regard to Cartesian coordinates x, y, the cable configuration follows 



amp; cos 0 (155) 

amp; sin 0. (156) 

The simultaneous integration of all four equations (T, 0, x, y) defines the cable 
configuration, and current dependency may be included, say U is a function 
of y. 

13.1.2 Critical Angle 

For very deep systems, the total weight of cable will generally exceed that the 
vehicle. This gives rise to a configuration in which the cable is straight for 
a majority of its length, but turns as necessary at the vehicle end, to meet 
the bottom boundary condition. In the straight part of the cable, normal 
weight and drag components are equalized. The uniform angle is called the 
critical angle 0c, and can be approximated easily. Let the relative importance 
of weight be given as 



dx 




ds 


amp] = 


dy 




ds 


amp; = 



Wn 

pCndU^' 
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Vohiclfti whidi are towe<I haw some similarities to the wliioles that have been 
discussed so far. For example, towed whicles are often streamlinnl, and usu- 
ally need good directional stal>ilily. Some towerl vehicles might have active 
lifting surfaces or llirusters for attitude control. Ou the other hand, if tliey 
are to ire supported Iry a cable, towerl vehicles may l>e quite heavy in water, 
and do not have to be self-propelled. The cable itself is an important factor in 
the behavior of the complete towerl system, aiul iu this swtiou. we concentrate 
ou cable mechanics more than vehicle characteristics, which can geuerally be 
handled with the same tools as other vehicles, i.e., slender-body theory, wing 
theory, linearizalian, ete.. Some basic guidelines for whicle design are given 
at the en<l of this section. 

Modern cables can easily exceed SOfiOm in length, ewn a heavy steel cable wdlh 
2cm diameter. The ealtles are geuerally circular hr croas section, and may 
carry power conductors and multiple communication channels (ltl>er optic). 
The extreme L/D ratio for these caldes oltviatea any bendhig stiffness effects. 
Calde systents come in a variety of configurations, and one main division may 
be marie .simply of the density of the calJe. Light-tether systems are character- 
Izeii by ueiitrally-lruoyant (or nearly so) cables, with either a minimaJ vehicle 
at the end. as in a towed array, or a vehicle capalde of maneuvering itself, such 
as a remotely-operated whicle. The towed array is a relatiwly high-velocity 
system that nominally streams itul horizontally behhnl the vessel. An ROV, 
on the other liarul, operates at low .speerl, and mast havr* large propiilsors to 
control the tether if there are currents. Heavy systems, in contrast, employ a 
heavy cable and possibly a heavy weight; the rationale is tliat gravity will tend 
to keep the calde vertical and make the deployment robust against currents 
and towing speed. Tlie heavy sy.stems will generally transmit surface motions 
and tensions to the towed whicJe much more easily than light-tether sy.stems. 
We will not discuss light systems specifically here, but rather look at heavy 
systems. Most of the analysis can be arlapteil to either case, howewr. 



13.1 Statics 

13.1.1 Force Balance 

For the purposes of deriving the static configuration of a cable in a flow, we 
assume for the moment that that it Is iuexlenslble. Tension anil hydrostatic 



SO that the condition d(j)/ds = 0 requires from the force balance 



5 cos (f)c 




0 . 



We are considering the case of 0 < 0c < Substituting sin^ (j)c = 1 — cos^ (f)c, 
we solve a quadratic equation and keep only the positive solution: 



cos 0c = a/5^ + 1 — 1. (157) 

In the case of a very heavy cable, S is large, and the linear approximation of 
the square root \/l + e ~ 1 + e/2 gives 



cos 0c amp] ~ 
0c amp; ~ 



amp: 



1 



71 

amp; - 



1 



(158) 



For a very light cable, 6 is small; the same approximation gives 
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cos 0c — 1 — ^ 0c — v^. 

The table below gives some results of the exact solution, and the approxima- 
tions. 



5 


amp; exact 


amp; S » 1 


amp; S « 1 


0.1 


amp; 0.44 


amp; - 


amp; 0.45 


0.2 


amp; 0.61 


amp; - 


amp; 0.63 


0.5 


amp; 0.91 


amp; 0.57 


amp; 1.00 


1.0 


amp; 1.14 


amp; 1.07 


amp; 1.41 


2.0 


amp; 1.33 


amp; 1.32 


amp; - 


5.0 


amp; 1.47 


amp; 1.47 


amp; - 



13.2 Linearized Dynamics 

13.2.1 Derivation 

The most direct procedure for deriving useful linear dynamic equations for a 
planar cable problem is to consider the total tension and angle as made up of 
static parts summed with dynamic parts: 



T{s,t)=f{s)+f{s,t) 

= 0(s) +0(s,t). 

We also write the axial deflection with respect to the static conhguration as 
p{s,t), and the lateral deflection q{s,t). It follows that 0 = dq/ds. Now 
augment the two static conhguration equations with inertial components: 



d‘^P dT df jjr . f7 7x 1 ^ 7 /rr 

m— amp; = amp; ^ ^ + 0) - \ ^ 



(m ^a)-^ amp; = amp; {T + T) \— + — j - Wn cos(0 + 0) 

^ rJn ■ ^ 

amp; amp; -pCnd I U sm 0 ~ j • 



+ 
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Here the material mass of the cable per unit length is m, and its transverse 
added mass is rua- Note that avoiding the drag law form n|n| again, we have 
implicitly assumed that U cos0 > \dp/dt\ and U sin0 > \dq/dt\. If it is not 
the case, say = 0, then equivalent linearization can be used for the quadratic 
drag. 

Now we perform the trigonometry substitutions in the weight terms, let 0 ~ 0 
for the calculation of drag, and substitute the constitutive (Hooke’s) law 



^ = EA^^. 
ds ds^ 

The static solution cancels out of both governing equations, and keeping only 
linear terms we obtain 



&^p 



[m 



m—w amp] = 
oH 

d‘^q 



d^p - dq - dp 

amp] EA^:— — Wn cos 4>— pCtdU cos 0— 

os^ os ot 

amp]f^ + EA^^ + fT^sin^^ - pCndU sin0^. 
os^ os os os ot 



The axial dynamics {p) couples with the lateral equation through the weight 
term —Wn cos 00. The lateral dynamics (g) couples with the axial through the 
term Td^/ds. An additional weight term WnSiiupdq/ds also appears. The 
uncoupled dynamics are both in the form of damped wave equations 



d‘^p , dp 

"'w + '"a 

dq 



d^q 

{m + ma)^ + h 



amp] = 



amp] = 



amp] EA 



d'^p 

ds'^ 



- d^q 

amp] T + Wn sin 0 
ds^ 



dq 

ds' 



where we made the substitution bt = pCtdU cos 0 and bn = pCndU sin 0. To a 
linear approximation, the out-of-plane vibrations of a cable are also governed 
by the second equation above. 

Because of light damping in the tangential direction, heavy cables easily trans- 
mit motions and tensions along their length, and can develop longitudinal 
resonant conditions (next section). In contrast, the lateral cable motions are 
heavily damped, such that disturbances only travel a few tens or hundreds of 
meters before they dissipate. The nature of the lateral response, in and out 
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of the towing plane, is a very slow, damped nonlinear filter. High-frequency 
vessel motions in the horizontal plane are completely missed by the vehicle, 
while low-frequency motions occur sluggishly, and only after a signihcant delay 
time. 





amp; axial 


amp; lateral 


wave speed 


j J? A 

amp; ^ 
amp; FAST 


amp; ^ 
amp; SLOW 


natural frequency 


amp; 


amp; 0 


mode shape 


amp; sine/cosine 


amp; Bessel function 


damping 


amp; Ct ^ 0(0.01) 


amp; 0„ ~ 0(1) 


disturbances 


amp; 


amp; 


travel down 


amp; YES 


amp; NO 


cable? 


amp; 


amp; 



13.2.2 Damped Axial Motion 

Mode Shape. The axial direction is of particular interest, since it is lightly 
damped and forced by the heaving of vessels in seas. Consider a long cable 
governed by the damped wave equation 



mp -I- btp = EAp" (159) 

We use over-dots to indicate time derivatives, and primes to indicate spatial 
derivatives. At the surface, we impose the motion 

p{L,t) = P cosuit, (160) 

while the towed vehicle, at the lower end, is an undamped mass responding to 
the local tension variations: 

EAdlpd ^ Mp(0,t). (161) 

OS 

These top and bottom behaviors comprise the boundary conditions for the 
wave equation. We let p{s,t) = p{s) cosuit, so that 
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~n I 

p + 



' moj"^ — ioobt 

'ea 



p = 0. 



(162) 



This admits the solution p(s) = ci cos ks + C2 sin ks, where 



k = 



I — iojE 



EA 



(163) 



Note that k is complex when 6^ 7^ 0. The top and bottom boundary conditions 
give, respectively, 



P amp] = amp] ci cos kL + C2 sin kL 
0 amp] = amp] Ci + hc2. 



where 5 = EAk/u'^M. These can be combined to give the solution 



In the case that M 
Pcos ks/ cos kL. 



S cos ks — sin ks 



p = P 

S cos kL 
0, the scalar S - 



00 , simplifying the result to p = 



Dynamic Tension. It is possible to compute the dynamic tension via T = 
EAp' . We obtain 



~ ^ 6 sinks + cos ks 

There are two dangerous situations: 

• The maximum tension is T + |T| and must be less than the working load 
of the cable. This is normally problematic at the top of the cable, where 
the static tension is highest. 

• If |T| > T, the cable will nnload completely and then reload with ex- 
tremely high impulsive forces. This is known as snap loading; it occurs 
primarily at the vehicle, where T is low. 



Natural Frequency. The natural frequency can be fonnd by letting bt = 0, 
and investigating the singularity of p, for which 6 cos kL = sin kL. In general. 
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kL « 1, but we find that a first-order approximation yields u = EA/ LM , 
which is only a correct answer if M » mL, i.e., the system is dominated by 
the vehicle mass. Some higher order terms need to be kept. We start with 
better approximations for sin() and cos(): 



5 





kLll- 



{kLf 



Employing the definition for 5, and recalling that = k'^EA/m, we arrive at 



ml / {kLf\ 

~W V 2 J 



{kiy 



{kLf 



If we match up to second order in /cL, then 



u = 

This has the familiar form of the square root of a stiffness divided by a mass: 
the stiffness of the cable is EA/L, and the mass that is oscillating is M+mL/2. 
In very deep water, the effects of mL/2 dominate; if pc = m/A is the density 
of the cable, we have the approximation 



EA/L 
M + mL/2 



1 





A few examples are given below for a steel cable with E = 200 x lO^Pa, and 
Pc = 7000/c5'/m^. The natural frequencies near wave excitation at the surface 
vessel must be taken into account in any design or deployment. Even if a 
cable can withstand the effects of resonance, it may be undesirable to expose 
the vehicle to these motions. Some solutions in use today are: stable vessels 
(e.g., SWATH), heave compensation through an active crane, a clump weight 
below which a light cable is employed, and an S-shaped length of cable at the 
bottom formed with flotation balls. 
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L = 500m 


amp; Un = I5.0rad/s 


1000m 


amp; 7.6rad/s 


2000m 


amp; 3.7rad/s 


5000m 


amp; 1.5rad/s 



13.3 Cable Strumming 

Cable strumming causes a host of problems, including obvious fatigue when 
the amplitudes and frequencies are high. The most noteworthy issue with 
towing is that the vibrations may cause the normal drag coefficient Cn to 
increase dramatically - from about 1.2 for a non-oscillating cable to as high 
as 3.5. This drag penalty decreases the critical angle of towing, so that larger 
lengths of cable are needed to reach a given depth, and the towed system lags 
further and further behind the surface vessel. The static tension will increase 
accordingly as well. 

Strumming of cables is caused by the proximity of a preferred vortex forma- 
tion frequency us to the natural frequency of the structure Un- This latter 
frequency can be obtained as a zero of the lightly-damped Bessel function 
solution of the lateral dynamics equation above. The preferred frequency of 
vortex formation is given by the empirical relation us = 27rSU/d, where S 
is the Strouhal number, about 0.16-0.20 for a large range of Re. Strumming 
of amplitude d/2 or greater can occur for 0.6 < oos/'^n < 2.0. The book by 
Blevins is a good general reference. 

13.4 Vehicle Design 

The physical layout of a towed vehicle is amenable to the analysis tools of 
self-propelled vehicles, with the main exceptions that the towpoint presents 
a large mean force as well as some disturbances, and that the vehicle can be 
quite heavy in water. Here are basic guidelines to be considered: 

1. The towpoint must be located above the vehicle center of in- water weight, 
for basic roll and pitch stability. 

2. The towpoint should be forward of the aerodynamic center, for towing 
stability reasons. 

3. The combined center of mass (material and added mass) should be longi- 
tudinally between the towpoint and the aerodynamic center, and nearer 
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the towpoint. This will ensure that high-frequency disturbances do not 
induce excessive pitching. 

4. The towpoint should be longitudinally forward of the center of in-air 
weight, so that the vehicle enters the water hns hrst, and self-stabilizes 
with U > 0. 

5. The center of buoyancy should be behind the in- water center of weight, 
so that the vehicle pitches downward at small U, and hence the net lift 
force is downward, away from the surface. 

Meeting all of these criteria simultaneously is no small feat, and the perfor- 
mance of the device is very sensitive to small perturbations in the geometry. 
For this reason, full-scale experiments are commonly used in the design pro- 
cess. 
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14 TRANSFER FUNCTIONS & STABILITY 

The reader is referred to Laplace Transforms in the section MATH FACTS for 
preliminary material on the Laplace transform. Partial fractions are presented 
here, in the context of control systems, as the fundamental link between pole 
locations and stability. 

14.1 Partial Fractions 

Solving linear time-invariant systems by the Laplace Transform method will 
generally create a signal containing the (factored) form 



F(.) = 



K[s ^i)(s -|- Z 2 ) ■ ■ ■ (<s -|- 
{s + Pi){s + P2) ■ ■ ■ {s + Pn) 



(166) 



Although for the moment we are discussing the signal T(s), later we will see 
that dynamic systems are described in the same format: in that case we call 
the impulse response G{s) a transfer function. A system transfer function is 
identical to its impulse response, since L{S(t)) = 1. 

The constants —Zi are called the zeros of the transfer function or signal, and 
—Pi are the poles. Viewed in the complex plane, it is clear that the magnitude 
of Y (s) will go to zero at the zeros, and to inhnity at the poles. 

Partial fraction expansions alter the form of Y (s) so that the simple transform 
pairs can be used to hnd the time-domain output signals. We must have 
m < n; if this is not the case, then we have to divide the numerator by the 
denominator as necessary to hnd a simple form. 



14.2 Partial Fractions: Unique Poles 

Under the condition m < n, it is a fact that Y (s) is equivalent to 



V(.) = 



Oi 



02 



s + pi 



P2 



+ 



S+Pn 



(167) 



in the special case that all of the poles are unique and real. The coefficient 
Qi is termed the residual associated with the Tth pole, and once all these are 
found it is a simple matter to go back to the transform table and look up the 
time-domain responses. 
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How to find aR A simple rule applies: multiply the right-hand sides of the 
two equations above by (s -|- pi), evaluate them at s = —pi, and solve for a*, 
the only one left. 



14.3 Example: Partial Fractions with Unique Real Poles 



G{s) 



s{s + 6) 



Since we have a pure delay and m = n, we can initially work with G{s)/se 
We have 



s -f 6 

(s + 4)(s- 1) 



Cli Qj2 

s -I- 4 s — 1’ 



giving 



(ii 



Cl2 



(s+6)(s+4) 

_(s+4)(s-l) 



s=— 4 



(s+6)(s-l) 

(s+4)(s-l) 



-I 5=1 



2 

5 



7 

5 



Thus 



L-\G{s)/se~^^) = + 

5 5 

9(i) = + + 

5 5 

The impulse response is needed to account for the step change at t = 2. Note 
that in this example, we were able to apply the derivative operator s after 
expanding the partial fractions. For cases where a second derivative must be 
taken, i.e., m > n + 1, special care should be used when accounting for the 
signal slope discontinuity at f = 0. The more traditional method, exemplihed 
by Ogata, may prove easier to work through. 

The case of repeated real roots may be handled elegantly, but this condition 
rarely occurs in applications. 
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14.4 Partial Fractions: Complex- Conjugate Poles 

A complex-conjugate pair of poles should be kept together, with the following 
procedure: employ the form 



F(.) 



hs + 62 

(s -Fpi)(s -Fps) 



+ 



03 



S+P3 



+ • • • , 



( 168 ) 



where pi = P 2 (complex conjugate). As before, multiply through by (s + 
Pi)(s + P 2 ), and then evaluate at s = —pi. 



14.5 Example: Partial Fractions with Complex Poles 



G(s) 



s + 1 



bis + 62 (U . 



's + 1' 

- ^ - s=-j 

1+J 

bi 

b2 

s + l 



[&1S + b2],^_j 

-bij + b2 — 
-1 

1; also 
03 = 1. 



Working out the inverse transforms from the table of pairs, we have simply 
(noting that C = 0) 



g(t) = — cost -|- sint -|- l(t). 

14.6 Stability in Linear Systems 

In linear systems, exponential stability occurs when all the real exponents of e 
are strictly negative. The signals decay within an exponential envelope. If one 
exponent is 0, the response never decays or grows in amplitude; this is called 
marginal stability. If at least one real exponent is positive, then one element 
of the response grows without bound, and the system is unstable. 
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14.7 Stability Poles in LHP 

In the context of partial fraction expansions, the relationship between stability 
and pole locations is especially clear. The unit step function 1 (t) has a pole at 
zero, the exponential has a pole at —a, and so on. All of the other pairs 
exhibit the same property: A system is stable if and only if all of the poles 
occur in the left half of the complex plane. Marginally stable parts correlate 
with a zero real part, and unstable parts to a positive real part. 

14.8 General Stability 

There are two dehnitions, which apply to systems with input u(t) and output 

y{t)- 

1. Exponential. If u(t) = 0 and y(0) = ijo, then \yi(t)\ < for hnite 

a and 7 > 0. The output asymptotically approaches zero, within a 
decaying exponential envelope. 

2. Bounded-Input Bounded-Output (BIBO). Ify(0) = 0, and \fi{t)\ < 
7,7 > 0 and hnite, then \yi(t)\ < a,a > 0 and hnite. 

In linear time-invariant systems, the two dehnitions are identical. Exponential 
stability is easy to check for linear systems, but for nonlinear systems, BIBO 
stability is usually easier to achieve. 
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15 CONTROL FUNDAMENTALS 

15.1 Introduction 

15.1.1 Plants, Inputs, and Outputs 

Controller design is about creating dynamic systems that behave in useful 
ways. Many target systems are physical; we employ controllers to steer ships, 
fly jets, position electric motors and hydraulic actuators, and distill alcohol. 
Controllers are also applied in macro-economics and many other important, 
non-physical systems. It is the fundamental concept of controller design that 
a set of input variables acts through a given “plant” to create an output. 
Feedback control then uses sensed plant outputs to apply corrective inputs: 



Plant 


Inputs 


Outputs 


Sensors 


Jet aircraft 


elevator, rudder, etc. 


altitude, hdg 


altimeter, GPS 


Marine vessel 


rudder angle 


heading 


gyrocompass 


Hydraulic robot 


valve position 


tip position 


joint angle 


U.S. economy 


fed interest rate, etc. 


prosperity 


inflation. Ml 


Nuclear reactor 


cooling, neutron flux 


power level 


temp., pressure 



15.1.2 The Need for Modeling 

Effective control system design usually benehts from an accurate model of 
the plant, although it must be noted that many industrial controllers can be 
tuned up satisfactorily with no knowledge of the plant. Ziegler and Nichols, for 
example, developed a general recipe which we detail later. In any event, plant 
models simply do not match real-world systems exactly; we can only hope to 
capture the basic components in the form of differential or integro-differential 
equations. 

Beyond prediction of plant behavior based on physics, the process of system 
identification generates a plant model from data. The process is often prob- 
lematic, however, since the measured response could be corrupted by sensor 
noise or physical disturbances in the system which cause it to behave in unpre- 
dictable ways. At some frequency high enough, most systems exhibit effects 
that are difficult to model or reproduce, and this is a limit to controller per- 
formance. 
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15.1.3 Nonlinear Control 

The bulk of this subject is taught using the tools of linear systems analysis. 
The main reason for this restriction is that nonlinear systems are difficult 
to model, difficult to design controllers for, and difficult overall! Within the 
paradigm of linear systems, there are many sets of powerful tools available. 
The reader interested in nonlinear control is referred to the book by Slotine 
and Li. 

15.2 Representing Linear Systems 

Except for the most heuristic methods of tuning up simple systems, control 
system design depends on a model of the plant. The transfer function de- 
scription of linear systems has already been described in the discussion of the 
Laplace transform. The state-space form is an entirely equivalent time-domain 
representation that makes a clean extension to systems with multiple inputs 
and multiple outputs, and opens the way to standard tools from linear algebra. 

15.2.1 Standard State-Space Form 

We write a linear system in a state-space form as follows 



X = Ax + Bu -I- Gw (169) 

y = Cx + Du + V 



where 



• a; is a state vector, with as many elements as there are orders in the 
governing differential equations. 

• is a matrix mapping x to its derivative; A captures the natural dy- 
namics of the system without external inputs. 

• i? is an input gain matrix for the control input u. 

• G is a gain matrix for unknown disturbance tc; w drives the state just 
like the control u. 

• y is the observation vector, comprised mainly of a linear combination of 
states Cx (where G is a matrix. 
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• Du is a direct map from input to output (usually zero for physical sys- 
tems). 

• u is an unknown sensor noise which corrupts the measurement. 




15.2.2 Converting a State-Space Model into a Transfer Function 

There are a number of canonical state-space forms available, which can create 
the same transfer function. In the case of no disturbances or noise, the transfer 
function (or transfer matrix) can be written as 

G(s) = 44 = C(sl - A)-^B + D, 
u{s) 

where I is the identity matrix with the same size as A 
holds for y{s)/w{s), and clearly y{s)/v{s) = I. 

15.2.3 Converting a Transfer Function into a State-Space Model 

It may be possible to write the corresponding differential equation along one 
row of the state vector, and then cascade derivatives. For example, consider 
the following system: 



(170) 

. A similar equation 



my'\t) + by\t) + ky{t) 



u' it) + u{t) (mass-spring-dashpot) 
s -1- 1 
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Setting X = [y\yY , we obtain the system 



dx 


—h/m 


—k/m 


X 


1/m 


dt 


1 


0 




0 



y = [ 11 ]:? 



Note specifically that dx 2 /dt = x\, leading to an entry of 1 in the off-diagonal 
of the second row in A. Entries in the C-matrix are easy to write in this case 
because of linearity; the system response to u' is the same as the derivative of 
the system response to u. 



15.3 PID Controllers 

The most common type of industrial controller is the proportional-integral- 
derivative (PID) design. If u is the output from the controller, and e is the 
error signal it receives, this control law has the form 



C(.) 



u{t) 

Ujs) 

E{s) 



kpe{t) + ki [ e{r)dT + kde{t), 
Jo 

ki 



kp h kd-s 



IH \-TdS 

TiS 



(171) 



where the last line is written using the conventions of one overall gain kp, plus 
a time characteristic to the integral part (r*) and and time characteristic to 
the derivative part (r^). 

In words, the proportional part of this control law will create a control action 
that scales linearly with the error - we often think of this as a spring-like action. 
The integrator is accumulating the error signal over time, and so the control 
action from this part will continue to grow as long as an error exists. Finally, 
the derivative action scales with the derivative of the error. The controller will 
retard motion toward zero error, which helps to reduce overshoot. 

The common variations are: P, PD, PI, PID. 



15.4 Example: PID Control 

Consider the case of a mass (m) sliding on a frictionless table. It has a per- 
fect thruster that generates force u(t), but is also subject to an unknown 
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disturbance d(t). If the linear position of the mass is y(t), and it is perfectly 
measured, we have the plant 



my" it) = u(t) + d{t). 

Suppose that the desired condition is simply y{t) = 0, with initial conditions 
1/(0) = yo and //'(O) = 0- 

15.4.1 Proportional Only 

A proportional controller alone invokes the control law u{t) = —kpy{t), so that 
the closed-loop dynamics follow 



my"{t) = -kpyit) + d{t). 

In the absence of d{t), we see that y{t) = //oCosy^f, a marginally stable 
response that is undesirable. 

15.4.2 Proportional-Derivative Only 

Let u{t) = —kpyit) — kdy'it), and it follows that 



my"{t) = -kpyit) - kdy'it) + d{t). 

The system now resembles a second-order mass-spring-dashpot system where 
kp plays the part of the spring, and kd the part of the dashpot. With an 
excessively large value for kd, the system would be overdamped and very slow 
to respond to any command. In most applications, a small amount of overshoot 
is employed because the response time is shorter. The kd value for critical 
damping in this example is 2^mkp, and so the rule is kd < 2^mkp. The 
result, easily found using the Laplace transform, is 



~^d 4 - 

y{t) = yoe^ 



kd 

sm ujdt 

2mud 



cos OJdt + 
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where Ud = ^Amkp — k'^/2m. This response is exponentially stable as desired. 
Note that if the mass had a very large amount of natural damping, a negative 
kd could be used to cancel some of its effect and speed up the system response. 
Now consider what happens if d{t) has a constant bias do- it balances exactly 
the proportional control part, eventually settling out at y(t = oo) = do/kp. To 
achieve good rejection of do with a PD controller, we would need to set kp very 
large. However, very large values of kp will also drive the resonant frequency 
Ud up, which is unacceptable. 

15.4.3 Proportional-Integral-Derivative 

Now let u{t) = —kpy{t) — ki jQy{r)dT — kdy'{t)\ we have 

my”{t) = -kpy{t) - ki [ y{T)dr - kdy(t) + d(t). 

Jo 

The control system has now created a third-order closed-loop response. If 
d(t) = do, a time derivative leads to 



my'”{t) -h kpy\t) + kiy{t) + kdy'\t) = 0, 
so that y{t = cx)) = 0, as desired, provided the roots are stable. 



15.5 Heuristic Tuning 

For many practical systems, tuning of a PID controller may proceed without 
any system model. This is especially pertinent for plants which are open-loop 
stable, and can be safely tested with varying controllers. One useful approach 
is due to Ziegler and Nichols (1942), which transforms the basic characteristics 
of a step response (e.g., the input is 1(f)) into a reasonable PID design. The 
idea is to approximate the response curve by a hrst-order lag (gain k and time 
constant r) and a pure delay T : 



G{s) ~ 



ke 



-Ts 



(172) 



rs -|- 1 

The following rules apply only if the plant contains no dominating, lightly- 
damped complex poles, and has no poles at the origin: 
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p 


kp - 


= l.Or/T 






PI 


kp - 


= 0.9r/T 


ki = f).27T/T^ 




PID 


kp - 


= 1.2r/T 


ki = 0.60r/T^ 


kd = 0.60r 



Note that if no pure time delay exists (T = 0), this recipe suggests the pro- 
portional gain can become arbitrarily high! Any characteristic other than a 
true hrst-order lag would therefore be expected to cause a measurable delay. 

15.6 Block Diagrams of Systems 

15.6.1 Fundamental Feedback Loop 

The topology of a feedback system can be represented graphically by consid- 
ering each dynamical system element to reside within a box, having an input 
line and an output line. For example, the plant used above (a simple mass) 
has transfer function P{s) = 1/ms^, which relates the input, force m(s), into 
the output, position y{s). In turn, the PD-controller has transfer function 
C(s) = kp + kds; its input is the error signal E{s) = —y{s), and its output is 
force u{s). The feedback loop in block diagram form is shown below. 




15.6.2 Block Diagrams: General Case 

The simple feedback system above is augmented in practice by three external 
inputs. The hrst is a process disturbance, which can be taken to act at the 
input of the physical plant, or at the output. In the former case, it is additive 
with the control action, and so has some physical meaning. In the second case, 
the disturbance has the same units as the plant output. 

Another external input is the reference command or setpoint, used to create a 
more general error signal e(s) = r{s) — y{s). Note that the feedback loop, in 
trying to force e(s) to zero, will necessarily make y{s) approximate r{s). 

The hnal input is sensor noise, which usually corrupts the feedback signal y{s), 
causing some error in the evaluation of e(s), and so on. Sensors with very poor 
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noise properties can ruin the performance of a control system, no matter how 
perfectly understood are the other components. 




15.6.3 Primary Transfer Functions 

Some algebra shows that 



e 1 

r 1 + PC 

y ^ PC 

r 1 + PC 

u _ C 

r 1 + CP 

e/r = S relates the reference input and noise to the error, and is known as the 
sensitivity function. We would generally like S to be small at low frequencies, 
so that the tracking error there is small, y/r = T is called the complementary 
sensitivity function. Note that S' + T = 1, implying that these two functions 
must always trade off; they cannot both be small or large at the same time. 
Other systems we encounter again later are the (forward) loop transfer function 
PC, the loop transfer function broken between C and P: CP, and 

e _ -P 

du 1 + PC 

]L = ^ 

du 1 + PC 

u _ -CP 
du 1 C P 



dy 1 + PC 



= S 
= T 
= U. 






108 



15 CONTROL FUNDAMENTALS 



]L 

dy 

U 

dy 

e 

n 

y 

n 

u 

n 



1 

l + PC 
-C 

l + CP 
-1 

l + PC 
-PC 
l + PC 
-C 

l + CP 



s 

-u 

-s 

-T 

-U. 



If the disturbance is taken at the plant output, then the three functions S', T, 
and U (control action) completely describe the system. This will in fact be 
the procedure when we address loopshaping. 




109 



16 MODAL ANALYSIS 

16.1 Introduction 

The evolution of states in a linear system occurs through independent modes, 
which can be driven by external inputs, and observed through plant output. 
This section provides the basis for modal analysis of systems. Throughout, we 
use the state-space description of a system with H = 0: 



X = Ax -f Bu 
y = Cx. 

16.2 Matrix Exponential 

16.2.1 Definition 

In the instance of an unforced response to initial conditions, consider the 
system 



X = Ax, x(t = 0) = y. 

In the scalar case, the response is x{t) = giving a decaying exponential 
if a < 0. The same notation holds for the case of a vector x, and matrix A: 



x{t) = where 

= 1 + At + + ■■■ 

is usually called the matrix exponential. 

16.2.2 Modal Canonical Form 

Introductory material on the eigenvalue problem and modal decomposition 
can be found in the MATH FACTS section. This modal decomposition of A 
leads to a very useful state-space representation. Namely, since A = VAV~^, 
a transformation of state variables can be made, x = Vz, leading to 
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z = A^ + V ^Bu (173) 

y = CVz. 

This is called the modal canonical form, since the states are simply the modal 
amplitudes. These states are uncoupled in A, but may be coupled through the 
input {V~^B) and output {CV) mappings. The modal form is numerically 
robust for computations. 

16.2.3 Modal Decomposition of Response 

Now we are ready to look at the matrix exponential in terms of its con- 
stituent modes. Employing the above form for A, we hnd that 



= 1 + At + 



2 ! 



+ ■■■ 



= V (^I + At+^-^ + --^W^ 
= 

n 

i=l 



In terms of the response to an initial condition y, we have 



i=l 

The product wfx is a scalar, the projection of the initial conditions onto the 
Tth mode. If x is perpendicular to Wi, then the product is zero and the i’th 
mode does not respond. Otherwise, the i’th mode does participate in the 
response. The projection of the i’th mode onto the states x is through the 
right eigenvector h). 

For stability of the system, the eigenvalues of A, that is, Aj, must have neg- 
ative real parts; they are in fact the poles of the equivalent transfer function 
description. 
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16.3 Forced Response and Controllability 

Now consider the system with an external input u: 



X = Ax + Bu, x(t = 0) = X- 

Taking the Laplace transform of the system, taking into account the initial 
condition for the derivative, we have 



sx{s) — X = Ax{s) + Bu{s) — 
x{s) = {si — A)~^x + ~ 

Thus {si — A)~^ can be recognized as the Laplace transform of the matrix 
exponential In the time domain, the second term then has the form of a 
convolution of the matrix exponential and the net input Bu\ 

x{t) = [ Bu{T)dT 

Jo 

n i-t 

= y^ Bu{r)dT. 

Suppose now that there are m inputs, such that B = [bi,b2, ■ ■ ■ ,bm]- Then 
some rearrangement will give 



n m 

x{t) = ^Vi^{wJbk) / 

i=i k=i 

The product wfbk, a scalar, represents the projection of the /c’th control chan- 
nel onto the i’th mode. We say that the dth mode is controllable from the 
/c’th input if the product is nonzero. If a given mode i has wjbk = 0 for all 
input channels k, then the mode is uncontrollable. 

In normal applications, controllability for the entire system is checked using 
the following test: Construct the so-called controllability matrix: 

M^=[B,AB,---,A^~^B]. (174) 

This matrix has size n x {nm), where m is the number of input channels. If 
Me has rank n, then the system is controllable, i.e., all modes are controllable. 
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16.4 Plant Output and Observability 

We now turn to a related question: can the complete state vector of the system 
be observed given only the output measurements y, and the known control u! 
The response due to the external input is easy to compute deterministically, 
through the convolution integral. Consider the part due to initial conditions 
X- We found above 

n 

i=l 

The observation is y = Cx (r channels of output), and writing 

r ->T 1 
Cl 

C = ■ . 

->T 

[ Cr J 

the /c’th channel of the output is 

n 

Vkit) = ^{c^Vi)e^^\wJx)- 

i=l 

The i’th mode is observable in the /c’th output if the product cj h) ^ 0. We 
say that a system is observable if every mode can be seen in at least one output 
channel. The usual test for system observability requires computation of the 
observability matrix: 

Mo = [C^, A^C^, • • • , {A^Y~^C^] . (175) 

This matrix has size n x (rn); the system is observable if Mo has rank n. 
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17 CONTROL SYSTEMS - LOOPSHAPING 

17.1 Introduction 

This section formalizes the notion of loopshaping for linear control system 
design. The loopshaping approach is inherently two-fold. First, we shape the 
open-loop transfer function (or matrix) P{s)C{s), to meet performance and 
robustness specihcations. Once this is done, then the compensator must be 
computed, from from knowing the nominal product P{s)C{s), and the nominal 
plant P{s). 

Most of the analysis here is given for single-input, single-output systems, but 
the link to multivariable control is not too difficult. In particular, absolute 
values of transfer functions are replaced with the maximum singular values of 
transfer matrices. Design based on singular values is the idea of L 2 -control, or 
LQG/LTR, to be presented in the next lectures. 



17.2 Roots of Stability — Nyquist Criterion 

We consider the SISO feedback system with reference trajectory r{s) and plant 
output y{s), as given previously. The tracking error signal is dehned as e(s) = 
r(s) —y{s), thus forming the negative feedback loop. The sensitivity function 
is written as 



r(s) l + P{s)C{sy 

where P{s) represents the plant transfer function, and C'(s) the compensator. 
The closed-loop characteristic equation, whose roots are the poles of the closed- 
loop system, is 1 -|- P{s)C{s) = 0, equivalent to P{s)C_{s) + P{s)C{s) = 
0, where the underline and overline denote the denominator and numerator, 
respectively. The Nyquist criterion allows us to assess the stability properties 
of a system based on P{s)C{s) only. This method for design involves plotting 
the complex loci of P{s)C{s) for the range u = [— cxd, cx)]. There is no explicit 
calculation of the closed-loop poles, and in this sense the design approach is 
quite different from the root-locus method (see Ogata). 
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17.2.1 Mapping Theorem 

We impose a reasonable assumption from the outset: The number of poles in 
P{s)C{s) exceeds the number of zeros. It is a reasonable constraint because 
otherwise the loop transfer function could pass signals with inhnitely high 
frequency. In the case of a PID controller (two zeros) and a second-order zero- 
less plant, this constraint can be easily met by adding a high-frequency rolloff 
to the compensator, the equivalent of low-pass hltering the error signal. 

Let F{s) = 1 -|- P{s)C{s). The heart of the Nyquist analysis is the mapping 
theorem, which answers the following question; How do paths in the s-plane 
map into paths in the F-plane? We limit ourselves to closed, clockwise{CW) 
paths in the s-plane, and the remarkable result of the mapping theorem is 
Every zero of F{s) enclosed in the s-plane generates exactly one CW encir- 
clement of the origin in the F{s) -plane. Conversely, every pole of F{s) enclosed 
in the s-plane generates exactly one CCW encirclement of the origin in the 
F{s)-plane. Since CW and CCW encirclements of the origin may cancel, the 
relation is often written Z — P = CW . 

The trick now is to make the trajectory in the s-plane enclose all unstable poles, 
i.e., the path encloses the entire right-half plane, moving up the imaginary 
axis, and then proceeding to the right at an arbitrarily large radius, back to 
the negative imaginary axis. 

Since the zeros of F{s) are in fact the poles of the closed- loop transfer function, 
e.g., S{s), stability requires that there are no zeros of F{s) in the right-half 
s-plane. This leads to a slightly shorter form of the above relation: 

p = CCW. (176) 

In words, stability requires that the number of unstable poles in F{s) is equal 
to the number of CCW encirclements of the origin, as s sweeps around the 
entire right-half s-plane. 

17.2.2 Nyquist Criterion 

The Nyquist criterion now follows from one translation. Namely, encirclements 
of the origin by F{s) are equivalent to encirclements of the point (—1 -|- Oj) 
by F{s) — 1, or P{s)C{s). Then the stability criterion can be cast in terms of 
the unstable poles of P(s)C(s), instead of those of F{s): 



P = CCW i — closed-loop stability 



(177) 
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This is in fact the complete Nyquist criterion for stability. It is a necessary and 
sufficient condition that the number of unstable poles in the loop transfer func- 
tion P{s)C{s) must be matched by an equal number of CCW encirclements 
of the critical point (— 1 -f Oj). 

There are several details to keep in mind when making Nyquist plots: 

• If neither the plant nor the controller have unstable modes, then the loci 
of P{s)C{s) must not encircle the critical point at all. 

• Because the path taken in the s-plane includes negative frequencies (i.e., 
the negative imaginary axis), the loci of P{s)C{s) occur as complex 
conjugates - the plot is symmetric about the real axis. 

• The requirement that the number of poles in P{s)C{s) exceeds the num- 
ber of zeros means that at high frequencies, P{s)C{s) always decays such 
that the loci go to the origin. 

• For the multivariable (MIMO) case, the procedure of looking at individ- 
ual Nyquist plots for each element of a transfer matrix is unreliable and 
outdated. Referring to the multivariable dehnition of S{s), we should 
count the encirclements for the function [det{I + P{s)C{s)) — 1] instead 
of P{s)C{s). The use of gain and phase margin in design is similar to 
the SISO case. 

17.2.3 Robustness on the Nyquist Plot 

The question of robustness in the presence of modelling errors is central to 
control system design. There are two natural measures of robustness for the 
Nyquist plot, each having a very clear graphical representation. The loci need 
to stay away from the critical point; how close the loci come to it can be 
expressed in terms of magnitude and angle. 

• When the angle of P{s)C{s) is —180°, the magnitude |P(s)C(s)| should 
not be near one. 

• When the magnitude |P(s)C(s)| = 1, its angle should not be —180°. 

These notions lead to dehnition of the gain margin kg and phase margin 7 for 
a design. As the hgure shows, the dehnition of kg is diherent for stable and 
unstable P{s)C{s). Rules of thumb are as follows. For a stable plant, kg >2 
and 7 > 30°; for an unstable plant, kg < 0.5 and 7 > 30°. 




116 



1 7 CONTROL SYSTEMS - LOOPSHAPING 





Stable P(s)C(s) 



Unstable P(s)C(s) 



17.3 Design for Nominal Performance 

Performance requirements of a feedback controller, using the nominal plant 
model, can be cast in terms of the Nyquist plot. We restrict the discussion to 
the scalar case in the following sections. 

Since the sensitivity function maps reference input r(s) to tracking error e(s), 
we know that |5'(s)| should be small at low frequencies. For example, if one- 
percent tracking is to be maintained for all frequencies below u = X, then 
|S'(s)| < 0.01, Vcu < A. This can be formalized by writing 

|Wi(s)^(s)| <1, (178) 

where W\{s) is a stable weighting function of frequency. To force S{s) to be 
small at low cu, Wi{s) should be large in the same range. The requirement 
|lFi(s)S'(s)| < 1 is equivalent to |lFi(s)| < |1 -|- P(s)C(s)|, and this latter 
condition can be interpreted as: The loci of P{s)C{s) must stay outside the 
disk of radius hFi(s), which is to be centered on the critical point (— 1 -|- Oj). 
The disk is to be quite large, possibly inhnitely large, at the lower frequencies. 

17.4 Design for Robustness 

It is a ubiquitous observation that models of plants degrade with increasing 
frequency. For example, the DC gain and slow, lightly-damped modes or zeros 
are easy to observe, but higher-frequency components in the response may 
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be hard to capture or even excite repeatably. Higher-frequency behavior may 
have more nonlinear properties as well. 

The effects of modeling uncertainty can be considered to enter the nominal 
feedback system as a disturbance at the plant output, dy. One of the most 
useful descriptions of model uncertainty is the multiplicative uncertainty: 

P{s) = (1 + A(s) 1 T 2 (s))P(s). (179) 

Here, P{s) represents the nominal plant model used in the design of the control 
loop, and P{s) is the actual, perturbed plant. The perturbation is of the 
multiplicative type, A(s)lT 2 (s)P(s), where A(s) is an unknown but stable 
function of frequency for which |A(s)| < 1. The weighting function HA("S) 
scales A(s) with frequency; 1 ^ 2 ( 5 ) should be growing with increasing frequency, 
since the uncertainty grows. However, HA(<s) should not grow any faster than 
necessary, since it will turn out to be at the cost of nominal performance. 

In the scalar case, the weight can be estimated as follows: since P/P — 1 = 
AHA, it will suffice to let \P/P — 1| < \W 2 \- 

Example: Let P = k/{s — 1), where k is in the range 2-5. We need to create 
a nominal model P = ko/{s — 1), with the smallest possible value of HA, which 
will not vary with frequency in this case. Two equations can be written using 
the above estimate, for the two extreme values of k, yielding ko = 7/2, and 
W 2 = 3/7. 

For constructing the Nyquist plot, we observe that 

P(.s)C(s) = (H- A(s)HA(s))P(s)C(s). The path of the perturbed plant could 
be anywhere on a disk of radius |HA(>s)P(s)C'(s)|, centered on the nominal loci 
P(s)C(s). The robustness condition is that this disk should not intersect the 
critical point. This can be written as 

|i + PC'! 

1 
1 

where T is the complementary sensitivity function. The last inequality is thus 
a condition for robust stability in the presence of multiplicative uncertainty 
parametrized with HA- 



> |W2PC'| < — ^ 

\W2PC\ 

^ |i + PC'!” 

> \W2Tl ( 180 ) 
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17.5 Robust Performance 

The condition for good performance with plant uncertainty is a combination 
of the above two conditions. Graphically, the disk at the critical point, with 
radius \Wi\, should not intersect the disk of radius 114^21^^1, centered on the 
nominal locus PC. This is met if 

llTi^l + \W2T\ < 1. (181) 

The robust performance requirement is related to the magnitude |T’C'| at dif- 
ferent frequencies, as follows: 

1. At low frequency, |hhiS'| ~ |fTi/PG|, since |PG| is large. This leads 
directly to the performance condition \PC\ > \Wi\ in this range. 

2. At high frequency, W 2 T\ ~ 1142^^1, since \PC\ is small. We must 
therefore have |PG| < l/|hh 2 |, for robustness. 

17.6 Implications of Bode’s Integral 

The loop transfer function PC cannot roll off too rapidly in the crossover 
region. The simple reason is that a steep slope induces a large phase loss, 
which in turn degrades the phase margin. To see this requires a short foray 
into Bode’s integral. For a transfer function H{s), the crucial relation is 

1 d 

angle{H{juJo)) = — —{ln\H{juj)) ■ ln{coth{\u\/2))du, (182) 

7T J—oo dy 

where v = ln{uj/uJo). The integral is hence taken over the log of a frequency 
normalized with ujq. It is not hard to see how the integral controls the angle: 
the function ln{coth{\v\/2)) is nonzero only near z/ = 0, implying that the 
angle depends only on the local slope d{ln\H\) / dv . Thus, if the slope is large, 
the angle is large. 

Example: Suppose H{s) = uJq/s^, i.e., it is a simple function with n poles at 
the origin, and no zeros; Uq is a hxed constant. It follows that \H\ = 
and ln\H\ = —nln{u/uo), so that d{ln\H\)/du = —n. Then we have just 

71 717T 

angle{H) = / ln{coth{\u\ / 2))du = — — . 

71 J-oo 2 



( 183 ) 
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This integral is trivial to look up or compute. Each pole at the origin clearly 
induces 90° of phase loss. In the general case, each pole not at the origin 
induces 90° of phase loss for frequencies above the pole. Each zero at the 
origin adds 90° phase lead, while zeros not at the origin at 90° of phase lead 
for frequencies above the zero. In the immediate neighborhood of these poles 
and zeros, the phase may vary signihcantly with frequency. 

The Nyquist loci are clearly susceptible to these variations is phase, and 
the phase margin can be easily lost if the slope of PC at crossover (where 
the magnitude is unity) is too steep. The slope can safely be hrst-order 
{—20dB /decade, equivalent to a single pole), and may be second-order 
{—AQdB / deeade) if an adequate phase angle can be maintained near crossover. 

17.7 The Recipe for Loopshaping 

In the above analysis, we have extensively described what the open loop trans- 
fer function PC should look like, to meet robustness and performance speci- 
hcations. We have said very little about how to get the compensator C, the 
critical component. For clarity, let the designed loop transfer function be re- 
named, L = PC. We will use concepts from optimal linear control for the 
MIMO case, but in the scalar case, it suffices to just pick 



C = L/P. (184) 

This extraordinarily simple step involves a plant inversion. 

The overall idea is to hrst shape L as a stable transfer function meeting the 
requirements of stability and robustness, and then divide through by the plant. 

• When the plant is stable and has stable zeros (minimum-phase), the 
division can be made directly. 

• One caveat for the well-behaved plant is that lightly-damped poles or 
zeros should not be cancelled verbatim by the compensator, because 
the closed-loop response will be sensitive to any slight change in the 
resonant frequency. The usual procedure is to widen the notch or pole 
in the compensator, through a higher damping ratio. 

• Non-minimum phase or unstable behavior in the plant can usually be 
handled by performing the loopshaping for the closest stable model, and 
then explicitly considering the effects of adding the unstable parts. In 
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the case of unstable zeros, we find that they impose an unavoidable 
frequency limit for the crossover. In general, the troublesome zeros must 
be faster than the closed-loop frequency response. 

In the case of unstable poles, the converse is true: The feedback system 
must be faster than the corresponding frequency of the unstable mode. 

When a control system involves multiple inputs and outputs, the ideas from 
scalar loopshaping can be adapted using the singular value. We list below 
some basic properties of the singular value decomposition, which is analogous 
to an eigenvector, or modal, analysis. Useful properties and relations for the 
singular value are found in the section MATH FACTS. 

The condition for MIMO robust performance can be written in many ways, 
including a direct extension of our scalar condition 

a{WiS) + a{W2T) < 1. (185) 

The open-loop transfer matrix L should be shaped accordingly. In the follow- 
ing sections, we use the properties of optimal state estimation and control to 
perform the plant inversion for MIMO systems. 
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18 LINEAR QUADRATIC REGULATOR 

18.1 Introduction 

The simple form of loopshaping in scalar systems does not extend directly 
to multivariable (MIMO) plants, which are characterized by transfer matrices 
instead of transfer functions. 

The notion of optimality is closely tied to MIMO control system design. Op- 
timal controllers, i.e., controllers that are the best possible, according to some 
hgure of merit, turn out to generate only stabilizing controllers for MIMO 
plants. In this sense, optimal control solutions provide an automated design 
procedure - we have only to decide what hgure of merit to use. The lin- 
ear quadratic regulator (LQR) is a well-known design technique that provides 
practical feedback gains. 

18.2 Full- St ate Feedback 

For the derivation of the linear quadratic regulator, we assume the plant to be 
written in state-space form x = Ax -|- Bu, and that all of the n states x are 
available for the controller. The feedback gain is a matrix K, implemented as 
u = —K{x — X desired)- The system dynamics are then written as: 

x= {A- BK)x H- KXdesired- (186) 

Xdesired represents the vector of desired states, and serves as the external input 
to the closed-loop system. The “A-matrix” of the closed loop system is {A — 
BK), and the “B-matrix” of the closed-loop system is K. The closed- loop 
system has exactly as many outputs as inputs: n. The column dimension 
of B equals the number of channels available in m, and must match the row 
dimension of K. Pole-placement is the process of placing the poles of [A — BK) 
in stable, suitably-damped locations in the complex plane. 

18.3 Dynamic Programming 

There are at least two conventional derivations for the LQR; we present here 
one based on dynamic programming, due to R. Bellman. The key observation 
is best given through a loose example: 

Suppose that we are driving from Point A to Point C, and we ask what is 
the shortest path in miles. If A and C represent Los Angeles and Boston, 
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for example, there are many paths to choose from! Assume that one way or 
another we have found the best path, and that a Point B lies along this path, 
say Las Vegas. Let X be an arbitrary point east of Las Vegas. If we were to 
now solve the optimization problem for getting from only Las Vegas to Boston, 
this same arbitrary point X would be along the new optimal path as well. 
The point is a subtle one: the optimization problem from Las Vegas to Boston 
is easier than that from Los Angeles to Boston, and the idea is to use this 
property backwards through time to evolve the optimal path, beginning in 
Boston. 

Example: Nodal Travel. We now add some structure to the above experi- 
ment. Consider now traveling from point A (Los Angeles) to Point D (Boston). 
Suppose there are only three places to cross the Rocky Mountains, Bi, B2, R3, 
and three places to cross the Mississippi River, Ci,C2,C3.^ By way of no- 
tation, we say that the path from A to B\ is AB\. Suppose that all of the 
paths (and distances) from A to the R-nodes are known, as are those from the 
R-nodes to the C-nodes, and the C-nodes to the terminal point D. There are 
nine unique paths from A to D. 

A brute-force approach sums up the total distance for all the possible paths, 
and picks the shortest one. In terms of computations, we could summarize that 
this method requires nine additions of three numbers, equivalent to eighteen 
additions of two numbers. The comparison of numbers is relatively cheap. 
The dynamic programming approach has two steps. First, from each R-node, 
pick the best path to D. There are three possible paths from Ri to D, for 
example, and nine paths total from the R-level to D. Store the best paths 
as BiD\opt, B2D\opt, B^Dlopt- This operation involves nine additions of two 

^Apologies to readers not familiar with American geography. 
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numbers. 

Second, compute the distance for each of the possible paths from A to D, 
constrained to the optimal paths from the B-nodes onward: ABi + BiD\opt, 
AB 2 + B 2 D\opt, or AB 3 + BsD\opt. The combined path with the shortest 
distance is the total solution; this second step involves three sums of two 
numbers, and total optimization is done in twelve additions of two numbers. 
Needless to say, this example gives only a mild advantage to the dynamic 
programming approach over brute force. The gap widens vastly, however, as 
one increases the dimensions of the solution space. In general, if there are s 
layers of nodes (e.g., rivers or mountain ranges), and each has width n (e.g., n 
river crossing points), the brute force approach will take {sn^) additions, while 
the dynamic programming procedure involves only (n^(s — 1) + n) additions. 
In the case of n = 5, s = 5, brute force requires 6250 additions; dynamic 
programming needs only 105. 

18.4 Dynamic Programming and Full-State Feedback 

We consider here the regulation problem, that is, of keeping Xdesired = 0. 
The closed-loop system thus is intended to reject disturbances and recover 
from initial conditions, but not necessarily follow y-trajectories. There are 
several necessary dehnitions. First we dehne an instantaneous penalty function 
l{x{t),u{t)), which is to be greater than zero for all nonzero x and u. The cost 
associated with this penalty, along an optimal trajectory, is 

COO 

J= l{x{t),u{t))dt, (187) 

Jo 

i.e., the integral over time of the instantaneous penalty. Finally, the optimal 
return is the cost of the optimal trajectory remaining after time t: 

COO 

V{x{t),u{t)) = J l{x{T),u{T))dr. (188) 

We have directly from the dynamic programming principle 

V {x{t),u{t)) = min {l{x{t),u{t))5t + V {x{t + 5t),u{t + 5t))} . (189) 

The minimization of V{x{t),u{t)) is made by considering all the possible con- 
trol inputs u in the time interval (t,t + St). As suggested by dynamic pro- 
gramming, the return at time t is constructed from the return at t -|- St, and 
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the differential component due to l{x,u). If V is smooth and has no explicit 
dependence on t, as written, then 



V {x(t + St),u(t + St)) 



V(x(t),u(t)) + — — — St + h.o.t. — (190) 
ox at 

dV 

V (x(t),u(t)) + —{Ax{t) + Bu{t))St. 



Now control input u in the interval (t,t + St) cannot affect V {x(t) , u(t)) , so 
inserting the above and making a cancellation gives 



f <91^ 1 

0 = min < l(x{t),u(t)) + -—(Axit) + Buit)) > . (191) 

“ [ ox J 

We next make the assumption that V (x, u) has the following form: 

V{x,u) = ^x'^Px, (192) 

where P is a symmetric matrix, and positive definiteff^ It follows that 



dV 

dx 

0 



x'^P — ^ 

min |/(a;, u) + x^P{Ax + Pm)| . 



(193) 



We finally specify the instantaneous penalty function. The LQR employs the 
special quadratic form 



l{x,u) = -x^Qx + -u^Ru, (194) 

where Q and R are both symmetric and positive definite. The matrices Q 
and R are to be set by the user, and represent the main “tuning knobs” for 
the LQR. Substitution of this form into the above equation, and setting the 
derivative with respect to u to zero gives 

^Positive definiteness means that x^Px > 0 for all nonzero x, and x'^ Px = 0 if a; = 0. 
®This suggested form for the optimal return can be confirmed after the optimal controller 
is derived. 
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0 = u^R + x^PB (195) 

u'^ = -x^PBR-^ 
u = —R~^B^Px. 

The gain matrix for the feedback control is thus K = R~^B^P. Inserting 
this solution back into equation 194, and eliminating u in favor of x, we have 

0 = ^x^Qx - ^x'^PBR-^B^P + x'^PAx. 

2^2 

All the matrices here are symmetric except for PA; since x'^PAx = x'^A'^Px, 
we can make its effect symmetric by letting 

x^PAx = ^x'^PAx + ^x'^A'^Px, 

leading to the hnal matrix-only result 

0 = Q + PA + A^P - PBR-^B^P. (196) 

This equation is the matrix algebraic Riccati equation (MARE), whose solution 
P is needed to compute the optimal feedback gain K. The MARE is easily 
solved by standard numerical tools in linear algebra. 

18.5 Properties and Use of the LQR 

Static Gain. The LQR generates a static gain matrix K, which is not a 
dynamical system. Hence, the order of the closed-loop system is the same as 
that of the plant. 

Robustness. The LQR achieves inhnite gain margin: kg = oo, implying that 
the loci of (PC) (scalar case) or {det{I + PC) — 1) (MIMO case) approach 
the origin along the imaginary axis. The LQR also guarantees phase margin 
7 > 60 degrees. This is in good agreement with the practical guidelines for 
control system design. 
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Output Variables. In many cases, it is not the states x which are to be 
minimized, but the output variables y. In this case, we set the weighting 
matrix Q = C^Q'C, since y = Cx, and the auxiliary matrix Q' weights the 
plant output. 

Behavior of Closed-Loop Poles: Expensive Control. When R >> 
C^Q'C, the cost function is dominated by the control effort m, and so the 
controller minimizes the control action itself. In the case of a completely stable 
plant, the gain will indeed go to zero, so that the closed-loop poles approach 
the open-loop plant poles in a manner consistent with the scalar root locus. 
The optimal control must always stabilize the closed-loop system, however, 
so there should be some account made for unstable plant poles. The expen- 
sive control solution puts stable closed-loop poles at the mirror images of the 
unstable plant poles. 

Behavior of Closed-Loop Poles: Cheap Control. When R« C^Q'C, 
the cost function is dominated by the output errors y, and there is no penalty 
for using large u. There are two groups of closed-loop poles. First, poles 
are placed at stable plant zeros, and at the mirror images of the unstable 
plant zeros. This part is akin to the high-gain limiting case of the root locus. 
The remaining poles assume a Butterworth pattern, whose radius increases to 
infinity as R becomes smaller and smaller. 

The Butterworth pattern refers to an arc in the stable left-half plane, as shown 
in the figure. The angular separation of n closed-loop poles on the arc is 
constant, and equal to 180° /n. An angle 90° /n separates the most lightly- 
damped poles from the imaginary axis. 
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19 KALMAN FILTER 



19 KALMAN FILTER 

19.1 Introduction 

In the previous section, we derived the linear quadratic regulator as an optimal 
solution for the full-state feedback control problem. The inherent assumption 
was that each state was known perfectly. In real applications, the measure- 
ments are subject to disturbances, and may not allow reconstruction of all the 
states. This state estimation is the task of a model-based estimator having 
the form: 



X = Ax + Bu + H{y — Cx) (197) 

The vector x represents the state estimate, whose evolution is governed by 
the nominal A and B matrices of the plant, and a correction term with the 
estimator gain matrix H. H operates on the estimation error mapped to 
the plant output y, since Cx = y. Given statistical properties of real plant 
disturbances and sensor noise, the Kalman Filter designs an optimal H. 




19.2 Problem Statement 

We consider the state-space plant model given by: 

X = Ax + Bu + Wi (198) 

y = Cx -f W 2 - 

There are n states, m inputs, and I outputs, so that A has dimension nxn, B 
is n X m, and G is / x n. The plant subject to two random input signals, Wi 
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and W 2 - W\ represents disturbances to the plant, since it drives x directly; 
W 2 denotes sensor noise, which corrupts the measurement y. 

An important assumption of the Kalman Filter is that Wi and W 2 are each 
vectors of unbiased, independent white noise, and that all the n + I channels 
are uncorrelated. Hence, if E{-) denotes the expected value. 



EiW,ih)Wi{t2f) 


— Vi^iti — 12 ) 


(199) 


E{W2ih)W2{t2f) 


= V25{E-t2) 


(200) 


E{Wi{t)W2{tf) 


~ ^nxl- 


(201) 



Here S{t) represents the impulse (or delta) function. Hi is an n x n diagonal 
matrix of intensities, and V 2 is an / x / diagonal matrix of intensities. 

The estimation error may be dehned as e = x — x. It can then be verihed that 



e = [Ax + Bu + W\] — [Ax + Bu + H{y — Cx)] (202) 

= {A- HC)e+{Wi- HW 2 ). 

The eigenvalues of the matrix A — HC thus determine the stability properties 
of the estimation error dynamics. The second term above, Wi + HW 2 is 
considered an external input. 

The Kalman hlter design provides H that minimizes the scalar cost function 

J = E{e^We), (203) 

where W is an unspecihed symmetric, positive definite weighting matrix. A 
related matrix, the symmetric error covariance, is dehned as 

E = E{ee^). (204) 

There are two main steps for deriving the optimal gain H. 

19.3 Step 1: An Equation for E 

The evolution of E follows from some algebra and the convolution form of e{t). 
We begin with 




130 



19 KALMAN FILTER 



S = E{ee^ + ee^) (205) 

= E[{A - HC)ee^ + {Wi - HW 2 )e^ + ee^{A^ - C^H^) + 
e{W[-W^H^)]. 

The last term above can be expanded, using the property that 

e{t) = - HW 2 {t)) dr. 

Jo 

We have 

E{e{W[ -W^H^)) = e^^-^^'^^E{e{0){W[ + 

Jo 

= - HW2{t)) 

{W[{t)-Wl{t)H^))dr 

= f - r) + HV 2 H^S{t - r)) dr 

Jo 

1 1 ^ 

= -Vi + -HV2H^. 

2 2 

To get from the hrst right-hand side to the second, we note that the initial 
condition e(0) is uncorrelated with — W^H^- The fact that Wi and HW 2 
are uncorrelated leads to the third line, and the hnal result follows from 

A 1 

S{t-T)dr = -, 

i.e., the written integral includes only half of the impulse. 

The hnal expression for E{e{W-[ — Wf^H'^)) is symmetric, and therefore ap- 
pears in Equation 205 twice, leading to 

t = {A- HC)T + - C^H'^) + El + HV2H^. (206) 

This equation governs propagation of the error covariance. It is independent 
of the initial condition e(0), and depends on the (as yet) unknown estimator 
gain matrix H . 
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19.4 Step 2: H as a Function of S 

We now make the connection between S = E{ee^) (a matrix) and J = 
E{e^We) (a scalar). The trace of a matrix is the sum of its diagonal ele- 
ments, and it can be verihed that 

J = E(e'^We) = trace{FW). (207) 

We now introduce an auxiliary cost function dehned as 

.J' = trace{T,W + AE), (208) 

where F is an n x n matrix of zeros, and A is an n x n matrix of unknown 
Lagrange multipliers. Note that since E is zero, J' = J, so minimizing J' solves 
the same problem. Lagrange multipliers provide an ingenious mechanism for 
drawing constraints into the optimization; the constraint we invoke is the 
evolution of S, Equation 206; 



J' = trace (fW + A(-S + AS - HCF + + W + HV2H^)) 

(209) 

If J' is an optimal cost, it follows that dJ' /dH = 0, i.e., the correct choice of 
H achieves an extremal value. We need the following lemmas, which give the 
derivative of a trace with respect to a constituent matrix: 

^trace{-AHCF) = 
uH 

^trace{-AFC'^ H^) = -AFC^ 
oH 

^trace{AHV2H^) = A^HV2 + AHV2. 

oH 

Proofs of the first two are given at the end of this section; the last lemma uses 
the chain rule, and the previous two lemmas. Next, we enforce A = A^, since 
the values are arbitrary. Then the condition dJ' / dH = 0 leads to 



= 2A(-SC'^ + HV 2 ), satished if 

= SC'^E2“^- 



0 

H 



( 210 ) 
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Duality of Linear Quadratic Regulator and Kalman Filter 



Linear Quadratic Regulator 


Kalman Filter 


X = Ax + Bu 


X = Ax -h Bu + Wi 




y = Cx + W 2 


u = —Kx 


X = Ax -f Bu + H{y — Cx) 


2J = J^{x'^Qx + vARu)dt 


J = E{e^'We) 


Q > 0,R > 0 


El P 0, V 2 ^ 0 


K = R-^B’^P 


H = 


PA + A^P + Q- PBR-^B^P = 0 


SA^ + AS + El - = 0 



Hence the estimator gain matrix H can be written as a function of S. Inserting 
this back into Equation 206, we obtain 

S = HS + TA^ + Hi - (211) 

Equations 210 and 211 represent the practical solution to the Kalman hltering 
problem, which minimizes the squared-norm of the estimation error. The 
evolution of S is always stable, and depends only on the constant matrices 
[A, C, Vi, V 2 ]. Notice also that the result is independent of the weighting matrix 
W, which might as well be the identity. 

19.5 Properties of the Solution 

The solution involves a matrix Riccati equation, like the LQR, suggesting a 
duality with the LQR problem. This is in fact the case, and the same analysis 
and numerical tools can be applied to both methodologies. 

The steady-state solution for S is valid for time-invariant systems, leading to 
a more common MARE form of Equation 211: 

0 = AS + SA^ + El - TC^V2~^CT. (212) 

The Kalman Filter is guaranteed to create a stable nominal dynamics A — HC, 
as long as the plant is fully state-observable. This is dual to the stability 
guarantee of the LQR loop, when the plant is state-controllable. Furthermore, 
like the LQR, the KF loop achieves 60° phase margin, and inhnite gain margin, 
for all the channels together or independently. 
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The qualitative dependence of the estimator gain H = on the other 

parameters can be easily seen. Recall that V\ is the intensity matrix of the 
plant disturbance, V2 is the intensity of the sensor noise, and S is the error 
covariance matrix. 

• A large uncertainty S creates large H, placing emphasis on the corrective 
action of the hlter. 

• A small disturbance Vi, and large sensor noise V2 creates a small H, 
weighting the model dynamics Ax + Bu more. 

• A large disturbance Vi, and small sensor noise V2 creates a large H, so 
that the hlter’s correction is dominant. 

The limiting closed-loop poles of the Kalman hlter are similar, and dual to 
those of the LQR; 



• V2 << Vi: good sensors, large disturbance, H » 1, dual to cheap- 
control problem. Some closed-loop poles go to the stable plant zeros, or 
the mirror image of unstable plant zeros. The remaining poles follow a 
Butterworth pattern whose radius increases with increasing V1/V2. 

• V2» Vp. poor sensors, small disturbance, H small, dual to expensive- 
control problem. Closed-loop poles go to the stable plant poles, and the 
mirror images of the unstable plant poles. 

19.6 Combination of LQR and KF 

An optimal output feedback controller is created through the use of a Kalman 
hlter coupled with an LQR full-state feedback gain. This combination is usu- 
ally known as the Linear Quadratic Gaussian design, or LQG. For the plant 
given as 



X = Ax + Bu -I- Wi 
y = Cx + W2, 



we put the Kalman Filter and controller gain G together as follows: 
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X = Ax + Bu + H{y — Cx) (213) 

u = —Kx. (214) 




d) = (sI-A')-^ 



There are two central points to this construction: 

1. Separation Principle: The eigenvalues of the nominal closed-loop sys- 
tem are made of up the eigenvalues of {A — HC) and the eigenvalues of 
{A — BK), separately. See proof below. 

2. Output Tracking: This compensator is a stand-alone system that, 
as written, tries to drive its input y to zero. It can be hooked up to 
receive tracking error e(s) = r(s) — y{s) as an input instead, so that it 
is not limited to the regulation problem alone. In this case, x no longer 
represents an estimated state, but rather an estimated state tracking 
error. We use the output error as a control input in the next section, on 
loopshaping via loop transfer recovery. 

19.7 Proofs of the Intermediate Results 

19.7.1 Proof that E{e^We) = trace(TW) 









d=i \t=i 

n n 
i=i j=i 



E{e^We) 
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the transpose of W being valid since it is symmetric. Now consider the diagonal 
elements of the prodnct HW: 



SW" = 



Siihhil + ^>12^21 + 



P‘21^12 + P‘22^22 + 



trace{TiW) = 

i=i j=i 

19.7.2 Proof that -^trace{—KHCTf) = — 



trace{AHB) = trace 



, the il'th element 



n I 

y~! ^ij HjkBki I 

Li=i fc=i 

n n I 

i=l j=l k=l 

where the second form is a snm over i of the ii’th elements. Now 



d 



dH 



jo^o 



d 



-trace{AHB) = Ajj^B, 



koi 



i=l 



— (BA)k^j^ 

= (BAf 



^—trace{AHB) 

oH 



= {BA) 
= A^B^ 



joko 

T 



□ 



19.7.3 Proof that j^trace{-AEC^ H^) = -AEC^ 

n 



trace{AH ^ ) = trace 



= trace 



,t=i 

n 

Y.AnH^ 



, the il'ih element 



i=l j=l 
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where the last form is a sum over the ii’th elements. It follows that 



d 



dK 



-trace{AH'^) = A 



'^ojo 



'^ojo 

-^—trace{AH'^) = A. 
oH 



□ 



19.7.4 Proof of the Separation Principle 

Without external inputs Wi and HA, the closed- loop system evolves according 
to 



d 

dt 




A -BK 

HC A-BK-HC 




Using the dehnition of estimation error e = x — x, we can write the above in 
another form; 




A- BK 
0 



BK 1 / a; 1 
A- HC \ \ e J ■ 



If A' represents this compound A-matrix, then its eigenvalues are the roots of 
det{sl — A') = 0. However, the determinant of an upper triangular block 
matrix is the product of the determinants of each block on the diagonal: 
det{sl — A') = det{sl — {A — BK))det{sI — {A — HC)) , and hence the separation 
of eigenvalues follows. 
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20 LOOP TRANSFER RECOVERY 

20.1 Introduction 

The Linear Quadratic Regulator(LQR) and Kalman Filter (KF) provide prac- 
tical solutions to the full-state feedback and state estimation problems, respec- 
tively. If the sensor noise and disturbance properties of the plant are indeed 
well-known, then an LQG design approach, that is, combining the LQR and 
KF into an output feedback compensator, may yield good results. The LQR 
tuning matrices Q and R would be picked heuristically to give a reasonable 
closed-loop response. 

There are two reasons to avoid this kind of direct LQG design procedure, 
however. First, although the LQR and KF each possess good robustness 
properties, there do exist plants for which there is no robustness guarantee 
for an LQG compensator. Even if one could steer clear of such pathological 
cases, a second problem is that this design technique has no clear equivalent 
in frequency space. It cannot be directly mapped to the intuitive ideas of 
loopshaping and the Nyquist plot, which are at the root of feedback control. 
We now reconsider just the feedback loop of the Kalman hlter. The KF has 
open-loop transfer function L{s) = C(p{s)H, where 0(s) = {si — A)~^ . This 
follows from the estimator evolution equation 



X = Ax -I- Bu + H{y — Cx) 

and the hgure. Note that we have not included the factor Bu as part of the 
hgure, since it does not affect the error dynamics of the hlter. 

As noted previously, the KF loop has good robustness properties, specihcally 
to perturbations at the output y, and further is amenable to output tracking. 
In short, the KF loop is an ideal candidate for a loopshaping design. Supposing 
that we have an estimator gain H which creates an attractive loop function 
L{s), we would like to hnd the compensator G(s) that establishes 



P{s)C{s) 

C(j){s)BC{s) 



^ C(f){s)H, or 
^ C(t){s)H. 



( 215 ) 



It will turn out that the LQR can be set up so that the an LQG-type com- 
pensator achieves exactly this result. The procedure is termed Loop Transfer 
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Recovery (LTR), and has two main parts. First, one carries out a KF design for 
H, so that the Kalman hlter loop itself has good performance and robustness 
properties. In this regard, the KF loop has sensitivity function S{s) = (/ + 
C(f){s)H)~^ and complementary sensitivity T(s) = (/ + C(f){s)H)~^C(f){s)H. 
The condition a(lFi(s)5'(s)) + a(lF 2 (-s)T(s)) < 1 is sufficient for robust per- 
formance with multiplicative plant uncertainty at the output. Secondly, we 
pick suitable parameters of the LQR design, so that the LQG compensator 
satishes the approximation of Equation 215. 

LTR is useful as a SISO control technique, but has a much larger role in 
multivariable control. 

20.2 A Special Property of the LQR Solution 

Letting Q = C^C and R = pi, where I is the identity matrix, we will show 
(roughly) that 



lim(VpiF) = WC, 

where K is the LQR gain matrix, and W is an orthonormal matrix, for which 
W^W = I. 

First recall the gain and Riccati equations for the LQR: 



K = R-^B^P 

0 = Q + PA + A'^P - PBR-^B^P. 

Now Q = C^C = C'^W'^WC = {WCYWC . The Riccati equation becomes 

0 = piWCfWC + pPA + pA^P - PBB^P = 0. 

In the limit as p — 0, it must be the case that P — 0 also, and so in this limit 
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piWCfWC 



wc 



^ pbb'^p 

= {B^PfB^P 
= {R-^B^PfRR{R-^B^P) 
= p^K^K — ^ 

y/pK. □ 



Note that another orthonormal matrix W could be used in separating K"’" 
from K in the last line. This matrix may be absorbed into W through a 
matrix inverse, however, and so does not need to be written. The result of the 
last line establishes that the plant must be square: the number of inputs (i.e., 
rows of K) is equal to the number of outputs (i.e., rows of C). 

Finally, we note that the above property is true only for LQR designs with 
minimum-phase plants, i.e., those with only stable zeros (Kwakernaak and 
Sivan). 



20.3 The Loop Transfer Recovery Result 

The theorem is stated as: If Xmip^Q^^K) = WC (the above result), with W 
an orthonormal matrix, then the limiting LQG controller C'(s) satisfies 



limP(s)G(s) = C(f)(s)H. 

The LTR method is limited by two conditions: 

• The plant has an equal number of inputs and outputs. 

• The design plant has no unstable zeros. The LTR method can be in fact 
be applied in the presence of unstable plant zeros, but the recovery is 
not to the Kalman hlter loop transfer function. Instead, the recovered 
function will exhibit reasonable limitations inherent to unstable zeros. 
See Athans for more details and references on this topic. 

The proof of the LTR result depends on some easy lemmas, given at the end 
of this section. First, we develop G(s), with the dehnitions 0(s) = {si — A)~^ 
and X(s) = (0“Hs) + HC)-^ = {sI-A + HC)~\ 
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C{s) = K{sl - A + BK + HC)-^H 

= K{X~^{s) + BK)~^H, then use Lemma 2 — 

= K{X{s) - X{s)B{I + KX{s)B)-^KX{s))H 

= KX{s)H - KX{s)B{I + KX{s)B)-^KX{s)H 

= (/ - KX{s)B{I + KX{s)B)-^)KX{s)H, then use Lemma 3 ^ 

= {I + KX{s)B)-^KX{s)H 
= {VpI + Vp^X{s)B)-^^KX{s)H. 

Next we invoke the result from the LQR design, with p — 0, to eliminate 
^K: 



limC'(s) = {WCX{s)B)-^WCX{s)H 

p^O 

= {CX{s)B)-^CX{s)H. 

In the last expression, we used the assumption that W is square and invertible, 
both properties of orthonormal matrices. Now we look at the product CX{s): 



CX{s) = C{SI-A + HC)-^ 

= C{(j)~^{s) + HC)~^, then use Lemma 2 — 

= C{(j){s) - (j)is)H{I + C(l){s)H)-^C(j){s)) 

= {I — C(j){s)H{I + C(j){s)H)~^)C(j){s), then use Lemma 3 
= (I + C(j){s)H)-^C(j){s). 

This result, reintroduced into the limiting compensator, gives 



limC'(s) = {{I + C(^{s)H)-^C(^{s)B)-\l + C(^{s)H)-^C4>{s)H 

p^O 

= {C(P{s)B)-^C(j){s)H 
= p-\s)C(j){s)H. 



Finally it follows that limp^o -P('S)C(s) = C(j){s)H, as desired. 
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20.4 Usage of the Loop Transfer Recovery 

The idea of LTR is to “recover” a Kalman filter loop transfer function L{s) = 

C(f){s)H, by using the limiting cheap-control LQR design, with Q = C^C and 

R = pi. The LQR design step is thus trivial. 

Some specihc techniques are useful. 

• Scale the plant outputs (and references), so that one unit of error in one 
channel is as undesirable as one unit of error in another channel. For 
example, in depth and pitch control of a large submarine, one meter of 
depth error cannot be compared directly with one radian of pitch error. 

• Scale the plant inputs in the same way. One Newton of propeller thrust 
cannot be compared with one radian of rudder angle. 

• Design for crossover frequency. The bandwidth of the controller is 
roughly equal to the frequency at which the (recovered) loop transfer 
function crosses over OdB. Often, the bandwidth of is a more intuitive 
design parameter than is, for example, the high-frequency multiplicative 
weighting IV 2 . Quantitative uncertainty models are usually at the cost 
of a lengthy identihcation effort. 

• Integrators should be part of the KF loop transfer function, if no steady- 
state error is to be allowed. Since the Kalman hlter loop has only as many 
poles as the plant, the plant input channels must be augmented with the 
necessary additional poles (at the origin). Then, once the KF design is 
completed, and the compensator C (s) is constructed, the integrators are 
moved from the plant over to the input side of the compensator. The 
tracking errors will accrue as desired. 

20.5 Three Lemmas 

Lemma 1: Matrix Inversion 



{A + BCD)-^ = R-i - A-^B{C~^ + D A~^ B)~^ D A~^ . 



Proof: 
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{A + BCD){A + BCD)-^ 
A{A + BCD)-^ 
{A + BCD)-^ 



I 

I - BCD{A + BCD)-^ 

^-1 _ A~^BCD{A + BCD)-^ 

- A~^BCD{I + A-^BCD)-^A~^ 

A-i - A-^B{D-^C~^ + A-^B)-^A~^ 

A~^ - A-^B{C~^ + DA-^B)-^DA~\ □ 



Lemma 2: Short Form of Lemma 1 



(X-i + BD)-^ = X - XB{I + DXB)-^DX 

Proof; substitute A = X~^ and C = I into Lemma 1. 

Lemma 3 



I-A{I + A)-^ = {I + A)~^ 



I -A{I + A)-^ 



{I + A){I + A)-^ - A{I + A)-^ 
{I + A - A){I + A)-^ 

(/ + A)-b □ 



Proof: 




143 



21 SYSTEM IDENTIFICATION 

21.1 Introduction 

Model-based controller design techniques, such as the LQR and LTR, require 
a plant model. The process of generating workable models from observed 
data is the goal of system identihcation. Good controllers usually have some 
reasonable robustness guarantees, which motivates identihcation with simple 
methods. We discuss in this section four fundamental but useful techniques 
for approximate system identihcation of single-input, single-output plants. It 
should be noted that the area of system identihcation is a very rich one, and 
that the methods are only a small subset of what is available. 

Except for the last approach, time-domain simulation, the methods are limited 
to linear models. If diherent inputs give diherent linear model coefficients, then 
it is likely that nonlinear terms are playing a role. The user then has the choice 
of ignoring the nonlinearity, for example if the operating point is controlled 
very closely, or developing a controller which takes specihc account of the it. 
In any event, simulations with the nonlinear plant should always be performed 
to assess the robustness of the control strategy. 

21.2 Visual Output froui a Siuiple luput 

For low-order plants which can tolerate impulse or step input, a great amount 
can be learned through step and impulse responses. The basic idea is to 
express what is observed as a time signal whose Laplace-domain equivalent 
can be recognized. As an example, consider the plant transfer function P{s) = 
k/{rs -fl), a first-order lag with gain k, and time constant r. We have y{s) = 
k/s{rs + 1) for the step input u{s) = 1/s, and therefore y{t) = k{l — 

The gain k is evident as the maximum value taken by the measured output. 
The time constant r is equal to the time required for y{t) to reach the value 
k{l — e~^) = 0.632A;. Similar estimates can be made for second-order systems, 
especially with the help of step functions parametrized on damping ratio (. 
Systems with order three or higher will usually be more difficult to assess with 
this visual technique. 

Example: Two raw step responses in heading were recorded for two dif- 
ferent vessels. The first vessel was strongly unstable and had to be powered 
down abruptly after one second; however, the smoothed data can still be £t- 
ted to the curve y(t) = 0.58(e* — 1). Give estimates of the two plant transfer 
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functions P{s), assuming that the step input was of magnitude one. 

The time trace of the hrst system looks like the Laplace transform pair: 



y{t) = - e-^') 

0 — a 



y{s) 



1 

(s + a) {s + h) 



for the values 6 = 0 and a = — 1. Thus for the experimental data, 



y{s) 



0.58 

s(s-l)‘ 



Since y{s) = P(s)u(s), where u{s) = 1/s, we have P{s) = 0.58/(s — 1). The 
second trace looks like a second-order response of the form 



y{s) ^ kul 

u{s) + 2(uJnS + up 

where gain k, undamped frequency u^ and ( are to be determined. First, 
we note that the steady value of y(t) is about 0.5, so let k = 0.5. Next, 
with respect to the steady value, the hrst overshoot is about 0.3, and the 
second overshoot (same side) is about 0.1. The ratio is often written as the 
logarithmic decrement 5 = /n(0. 3/0.1) ~ 1.10, so that the damping ratio is 
simply = 5/2ti ~ 0.175. Finally, the damped natural period is about 1.7s, 
leading to the damped natural frequency Ud = 27t/ 1.7 ~ 3.70rad/s. The 
undamped natural frequency Un is related to Ud as follows: 
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UJn 



(^d 



~ 3.75rad/s. 



Inserting these into the template above, we have 



+ 1.28s + 14.r 

21.3 Transfer Function Estimation — Sinusoidal Input 

The main idea of transfer function estimation is that P{s) can be estimated by 
simply dividing the measured output by the measured input, in the frequency 
domain: 



P(s) = ^ (216) 

m(s) 

In applications, the input signal u{s) is known quite accurately because it is 
generated by a computer. Two sources of error can corrupt the output signal 
y{s), however: real disturbances and sensor noise. In some cases, disturbances 
may be band-limited (e.g., water waves), and if these occur in a frequency 
range far away from the dominant dynamics, the estimated transfer function 
approach will succeed. A similar argument holds for sensor noise, which in 
many devices is negligible for low frequencies. A sensor which has noise in the 
frequency range of the system dynamics is problematic for obvious reasons. 
This section discusses the use of periodic inputs to create a Bode plot, while 
the next section generalizes to broadband input. Bode plots are hgures of 
transfer function magnitude and phase as functions of frequency, which can 
be either parameterized in terms of poles, zeros, and a gain, or used directly 
in a loopshaping or Nyquist plot approach. 

Certain plants which cannot admit a step or impulsive input will tolerate a 
sinusoidal input. The idea then is to drive the plant with a narrow-band, i.e., 
periodic, input signal u{s). For each such test at a specific frequency, compute 
the magnitude and phase relating the input to the output. Conduct as many 
tests as are necessary to build a Bode plot of the plant estimate P{s). 

Example: Sinusoidal rudder angles trajectories, of lOdeg amplitude (|M(t)| = 
lOdeg), were implemented on a vessel operating near its cruising speed. The 
data are plotted, along with the magnitude and phase of the transfer function 
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p _ 0.093 

“ s2 + 0.83s + 0.21’ 

which holds reasonably well at low frequencies. The phase angle of y{s) is taken 
with respect to the input signal u{s). Note that the experimental magnitude 
and phase in this example deteriorate at the higher frequencies. This is a 
property of almost all physical systems, and an indicator that the plant model 
cannot be trusted above a certain frequency range. 




21.4 Transfer Function Estimation — Broadband Input 

Many times one needs to deal with experimental data that is broadband, either 
as the result of a closed-loop run, or of significant disturbances. In this case, 
frequency-domain analysis can still be used, but elements of spectral analysis 
are necessary. 
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21.4.1 Fourier Transform of Sampled Data 

For the purpose of analyzing transfer functions, the discrete Fourier transform 
(DFT) is a standard, optimized tool. It operates on a data vector x{n) of 
length N: 



N-l 

X{m) = Y. (217) 

n=0 

for m = [1, A^]. We review several points for dealing with the DFT calculation: 

• The m’th DFT point is a summation of complex unit-magnitude vectors 

times the original data {x{n)). Dividing the DFT result by 
N returns the signal to its real magnitude. 

• The DFT generates a vector of N complex points as outputs. These 
correspond to the frequency range 



27rm — 1 
^ hi N 



( 218 ) 



i.e., the first frequency is 0, and the last frequency is slightly less than 
the sampling rate. 



• The sampling theorem limits our usable frequency range to only one-half 
of the sampling rate, called the Nyquist rate. The points higher than 
Ti /dt correspond to negative frequencies, and a complex-conjugacy holds: 
X(2) = X*(m),X(3) = X*(m — 1), and so on. What happens near the 
Nyquist rate depends on whether N is odd or even. 

An additional multiplication of the DFT result by 2 will give peaks which 
are of about the right magnitude to be compared with the time-domain 
signal. With this scaling, the value W(l) is twice the true DC value. 

• The Nyquist rate depends only on the sampling time step dt, but the 
frequency vector accompanying the DFT can be made arbitrarily long 
by increasing the number of data points. To improve the frequency 
resolution of the DFT, a common approach is to zero-pad the end of the 
real data. 
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• The data x{k) should be multiplied by a smoothing window w{k), for 
two reasons. First, if the window goes to a small value near its tails, 
Gibb’s effect (the spectral signature of a discontinuity between a;(l) and 
x{N)) can be minimized. The second rationale for windowing is based 
on the discrete nature of the transform. Because the DFT provides 
frequency data at only N/2 unique frequencies, there is the possibility 
that a component in the real data lies in between one of these DFT 
frequencies. The DFT magnitude at a specihc u{m) is in fact an average 
of continuous frequencies from the neighborhood of the point. Smoothing 
windows are generally chosen to achieve a tradeoff of two conflicting 
properties relating to the average: the width of the primary lobe (wide 
primary lobes bring in frequencies that belong to the adjacent bins u{m— 
1) and + 1)), and the magnitude of the side lobes (tall side lobes 
bring in frequencies that are far away from u{m)). 

There are many smoothing windows to choose from; no windowing at 
all is usually referred to as applying a rectangular window. One of the 
simplest is the Hann, or cosine, window: 

1 

w{k) = - 

• The DFT of a given signal is subject to bias and variance. These can 
be reduced by taking separate DFT’s of sub-sections of data (perhaps 
zero-padded to maintain frequency resolution), and then averaging them. 
It is common to use segments which do not overlap, but it is not a 
necessity. This averaging idea is especially important in transfer function 
verihcation; see below. 

• If possible, the DFT should be run on a number of samples N which is 
a multiple of a large power of two. If N is prime, the DFT will take a 
long time to execute. 

The above steps, all of which are available through the Matlab function 
spectrum 0 for example, will ensure a fair spectral analysis of a time series. 

21.4.2 Estimating the Transfer Function 

In continuous time, when the plant input y(t) and output u(t) are transformed 
into frequency space, the estimated transfer function follows from 
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P(s) = ^ (219) 

m(s) 

As noted above, spectral analysis of a signal benefits by using multiple seg- 
ments and averaging; we now want to include the same approach in our esti- 
mation of P{s). The procedure is: 

• For the p’th segment of data, compute the transfer function estimate 



_ Yp{m)Y*{m) 
Up{m)U*{m) ' 



• For the p’th segment of data, compute two covariances and a cross- 
covariance: 



rL(m) = U,(m)U;(m) 

KM) = Fp(m)y;(m) 

KM) = Y,(m)u;(m). 

• Construct average values of the transfer function estimate and the other 
quantities: 



P(m) = avgp{Pp{m)) 

f'uuim) = avgp(Tl^{m)) 
Tyy{m) = avgp{Tly{m)) 
Tyu{m) = avgp(TPy^{m)) 



• Compute the coherence function, which assesses the quality of the hnal 
estimate P{m)\ 



^ yu{m)V*y^{m) 
V yy{m)T uu{m)' 



Coh{m) 
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If the coherence is near zero, then the segmental cross-covariances 
are sporadic and have cancelled out; there is no clear relation between 
the input and the output. This result could be caused by either dis- 
turbances or sensor noise, both of which are reasonably assumed to be 
random processes, and uncoupled to the input signal. Alternatively, if 
the coherence is near one, then the cross-covariances are in agreement 
and a real input-output relationship exists. With real data, the coher- 
ence will deteriorate at high frequencies and also at any frequency where 
disturbances or noise occur. 

21.5 Time-Domain Simulation 

The time-domain simulation approach tweaks the parameters of a simulation 
so that its output matches the observed output. The method has its main 
strength in the fact that it applies to any model that can be simulated, in- 
cluding those of high order and with signihcant nonlinearities. On the other 
hand, the method is computationally expensive and gives no guarantee of a 
useful solution, or even of convergence. 

At the outset, we need to come up with some structure of the plant model. 
This can be based on physics in many cases. Consider, for example, the case 
of a mass mounted on a spring and a dashpot, driven by the input force u(t). 
A fair guess for the actual dynamics has the form my” + by' + ky = u{t), and 
we plan to look through the three-dimensional parameter space 9 = [m, b, k]. 
The simulation operation can be written this way: y(t) = G{9,u(t)), and the 
system identihcation problem is to minimize \ \y(t) — yobs(t)\\, say, where || ■ || 
here indicates the Euclidean norm. For a given parameter vector 6, running the 
simulation generates a new y, and computing the norm gives a scalar measure 
of goodness. 

Since the normed error is a complicated function of both u{t) and 6, the mini- 
mization must proceed iteratively. The Nelder-Meade simplex method is easy 
to use, and can be invoked with the Matlab function fmins. As an example, 
the three programs listed comprise a working Matlab set for identihcation of 
a hrst-order, nonlinear system. Some notes on use: 

• In this example, the same simulation generates the “observed” data and 
the simulated response. Since the program simulate always uses the 
global variable theta as the parameters, we must be careful about setting 
theta in the calling programs. 
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• After a simulation run is complete, the data is interpolated to the same 
time scale as the observed data, in order to compute the error. 

• The initial guess for theta is a random vector; the Simplex method will 
take over from this point. In many instances, however, theta is roughly 
known, and a better starting value can be given. 

• The Simplex method may head into invalid parameter space, e.g., nega- 
tive mass. The error calculation, however, can be easily augmented by 
a term which penalizes invalid parameters, e.g., 

err = err + 1000* (1-sign (mass) ) . 

• There is no guarantee that a global minimum will be found or even exists. 
Starting from different initial guesses for 0 may help hnd better results, 
but we are still at the mercy of the minimization algorithm, and a very 
complicated function. 



0 /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 



clear all ; clear global ; 
global u_obs t_obs y_obs dt theta ; 

dt = .3 ; % time step 

theta =[12] ; % true parameter vector 
t_obs = 0:dt:20*dt ; "Z, observed time vector 
u_obs = ones (length(t_obs) , 1) ; % observed input 
[t_raw, y_raw] = ode45 (^ simulate ^ , [0 max(t_obs)] , 0) ; 
y_obs = spline (t_raw, y_raw, t_obs) ; % observed output 

[theta_f inal] = fmins('get_err' , randn(2,l)) ; 

disp(sprintf ( 'final theta(l) : 7og.', theta_f inal (1) ) ) ; 

disp(sprintf ( 'final theta(2) : 7g.', theta_f inal (2) ) ) ; 

theta = theta_final ; 

[t_raw, y_raw] = ode45 (' simulate ' , [0 max(t_obs)] , 0) ; 
y_sim = spline (t_raw, y_raw, t_obs) ; 
figure(l) ; clf ; hold off ; 

plot(t_obs, y_obs, t_obs, y_sim, t_obs, u_obs) ; 
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0 / 0 / 0 / 0 / 0 / 0 / 



0 /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 



function [err] = get_err (theta_arg) ; 
global y_obs t_obs theta ; 

theta = theta_arg ; 

[t_raw, y_raw] = ode45 simulate ^ , [0 max(t_obs)] , 0) ; 
y_sim = spline (t_raw, y_raw, t_obs) ; 
err = norm(y_sim - y_obs) ; 
dispCsprintf (' error : %f . ’ , err)) ; 



0 / 0 / 0 / 0 / 0 / 0 / 



0 /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 



function [ydot] = simulate (t,y) ; 

global u_obs dt theta ; 

k = theta(l) ; 
tau = theta(2) ; 

ind = floor (t/dt) + 1 ; % we have to choose which u_obs to use: 
7o a zero-order hold as implemented. 

ydot = ( k*u_obs(ind) - y'S ) / tau ; 



0 / 0 / 0 / 0 / 0 / 0 . 



0 /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 
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The bulk of our discussion on maneuvering and control has assumed that the 
necessary system states can be measured. The marine engineer is in fact faced 
with choices between many different basic sensor packages, notably compasses, 
paddle wheels, inertial navigation units, rate gyros, and depth guages, for ex- 
ample. These listed sensors are self-contained and rely primarily on the phys- 
ical properties of the natural environment. There is also a class of distributed 
sensor systems; these generally involve an array of communicating elements, lo- 
cated remotely from the vehicle. We present the fundamental concepts behind 
two methodologies in this second class: the global positioning system (GPS) 
and acoustic navigation, both of which can provide high-accuracy absolute 
Cartesian navigation. 

22.1 Acoustic Navigation 

Consider a transponder A, which can transmit an acoustic signal, and also 
measure, with microsecond accuracy, the time to receive a reply. Next, place 
a responder B at a distance R away from A; the job of B is just to transmit 
a signal whenever it receives one, with a (short) predictable response time R. 
Thus, the elapsed time T between a tranmission and consequent reception at 
A is 



T — Tt ‘2.R/ c^, 

where c^, is the speed of sound in water, about 1450m/s. The range R follows 
by inversion: 



^ — ■ 

Suppose that the location of B is known; then a given measurement of R 
places A on a sphere around B. It is a case of three unknowns (x, y, z) and one 
equation: 



{xa - xbY + {va - VbY + {za - ZbY = R^. 




154 



22 CARTESIAN NAVIGATION 



The introduction of a another responder C (typically listening and responding 
at a different frequency than B) places A on the intersection of two spheres, 
i.e., a circle. There are two equations, but still three unknowns. When A, 
B, and C he in a nearly horizontal plane, then the intersection circle lies in a 
vertical plane; the addition of a depth sensor to our suite would allow us to pin 
A’s location at one of two points on the circle. Finally, if we know which side A 
is on, and do not allow for abrupt crossovers, then we have a functional set of 
measurements for acoustic navigation with just two responders. The baseline is 
the line connecting responders B and C; when A is near this baseline, positional 
accuracy will be very poor since the two solution spheres are tangent. 

Better and better performance can be obtained by increasing the number of 
responders, and consequently of the baselines and spheres. With three re- 
sponders, for example, the intersection of a sphere (responder D) and a circle 
(responders B and C) is two points. Here there are three equations and three 
unknowns, but the nonlinearity of the equations leads to the non-uniqueness 
in the solution. A fourth responder or a depth transducer would be needed to 
completely constrain the solution. 

The above discussion is a minimum conceptual explanation of acoustic naviga- 
tion. There are many other pieces to the approach, including an account of the 
variation of sound speed with water depth, obtaining the Cartesian loca- 
tions of the responder network, and handling various geometric conhgurations 
that give rise to poor or degenerate solutions. 

There are two common conhgurations used for acoustic navigation, named for 
the length of the baselines relative to the target (transponder A) range. 

(Ultra) Short-Baseline The geometry of SBL or USBL puts very short 
baselines between the responders compared to the target distance. For in- 
stance, the hxed net is often attached to a vessel or other structure, with base- 
line lengths on the order of lO-lOOm. Typical frequencies in use are around 
lOOkHz, with a working range of 100-500m to the target. The wavelength of a 
lOOkHz signal is about 1cm and 5cm is a reasonable estimate of the accuracy 
for these systems. 

Long-Baseline Long-baseline systems typically involve a larger responder 
net, and the target ranges are similar to the baseline lengths. Very large 
systems utilize frequencies of 10-15kHz (for a placement accuracy around 2- 
5m), and may have ten-kilometer baselines. Frequent sources of error in long- 
baseline systems are variations of c^, and also in multipath. In the latter 
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Figure 6: Left: General configuration of a long-baseline acoustic system, with 
target transponder A and fixed responders B and C. The position Ai is on 
the baseline and has poor accuracy in the direction normal to the baseline; in 
contrast, A3 is well-posed. A3 is an additional solution to the two-responder 
problem. Right: General conhguration of a short-baseline system. 



condition, false signals can be caused by reflections off the seafloor and the 
surface; these signals are typically eliminated by rejecting those receptions 
which are outside a very tight and slowly moving window on travel time T. 
Sometimes, bottom topography can shield a direct acoustic path, and the only 
receptions available are via multipath! 

22.2 Global Positioning System (GPS) 

The GPS system is the most powerful system publicly available for absolute 
position reference above water. It is similar in concept to the acoustic naviga- 
tion systems described above, bnt with one fnndamental difference: only the 
one-way travel time from each satellite is measnred at the target. The accu- 
racies achieved with GPS are therefore strongly dependent on the accnracy of 
time-keeping. 
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The core of the system is a network of 24 satellites arranged in six planes (4 
per plane) inclined at 55 degrees from the equatorial plane, at an altitude of 
about 20,000/cm. The altitude is chosen carefully to correspond with a 12-hour 
orbit; geosynchronous orbit is much higher, at around 40,000/cm. The orbit 
lines over a non-rotating earth are thus near-sinusoids that reach lattitudes 
of 55°N and 55°S, and the six lines are spaced 60° apart in longitude. With 
reference to a rotating Earth, the spatial frequency is doubled, so that the 
longitudinal distance between tracks is 30 degrees, or 1800 nautical miles. 
Each satellite transmits a regular signal which contains, among other items, the 
time of its transmission, and the three-dimensional location of the satellite (the 
ephemeris). It is from these two pieces of information, from many satellites, 
that triangulation and navigation are performed at the target. 

For the ideal case that the time bases of the target and the satellites are 
exactly aligned, triangulation of the type described for acoustic navigation is 
possible. In practical terms, the speed of light c = 3 x llRm/s implies that 
a Ins timing error will cause a 0.3m range error. For this reason, each of 
the satellites carries four atomic clocks on board, good to Ins per day. The 
time base is continually monitored and updated from a ground station, and 
accurate to within Ins. The ephemeris of each satellite is also monitored and 
updated from the ground station, using a combination of least-squares analysis 
of past data (1 week), and a Kalman hlter to predict the future ephemeris. 
The accuracy of satellite position is better than 10m. 

We now come to the last thorny issue: how to make target receivers that don’t 
require atomic clocks! The solution is very clever, and illustrated for the case 
of two-dimensions; the arguments for three dimensions are the same. In two 
dimensions, the range measurements from two sources locate the target on one 
of two points (the intersection of two circles). Suppose that the target clock 
is too slow by an amount t, so that the range estimate is 



r f = c{t + t). 

The estimated range circles are too large, and the estimated location of the 
target too far from the baseline. Introduce a third satellite now so that the 
intersections of the circle pairs now occur at six points. Three of these are 
very close, and indicate the approximate true solution. The trick is to hud the 
correction for t that puts the three close solutions onto a single point. This 
single point is near the centroid of the three approximate points, and represents 
the best position solution. The target time base can be kept up to date by 
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performing this check for every set of signals. Hence, the two-dimensional 
timing problem is solved with three satellites; the three-dimensional timing 
problem is solved with four satellites. 

The position specihcation for GPS is 25m, at the 95’th percentile. This is a re- 
markable feat, given that it is an absolute measure over the entire surface and 
atmosphere of our 6000/cm-radius Earth. Major sources of error include: clock 
base and satellite navigation 2.5m), ionosphere and troposphere electro- 
magnetic variations (~ 2.5m), and inaccuracies in the receiver and multipath 
2m). GPS is also subject to selective availability or S/A, the addition of 
a slow-varying random component in the satellite ephemeris data. Selective 
availability is controlled by the United States Department of Defense, and de- 
grades the position specihcation to 100m (typical). The advent of differential 
GPS solved most of the S/A concerns of American allies, by providing high- 
accuracy navigation in local areas. The idea here is that a stationary target 
can detect the effects of S'/A (since it is not moving) and then transmit the 
corrections over a small geographical area. Many current GPS receivers are 
capable of decoding these local corrections, which are then applied to their 
own satellite navigation processing. Differential GPS typically provides l-2m 
accuracy. 

Finally, as with the acoustic navigation systems described, an independent 
altitude measurement will enhance the accuracy of GPS, essentially reducing 
a three-dimensional to a two-dimensional problem. 
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24 PROBLEMS 

Many of the following problems use data presented by Fossen and in Principles 
of Naval Architecture. 

1. Inertial dynamics. The gas generator turbine of a LM-2500 marine 
engine rotates at p = 10, OOOrpm (revolutions per minute), and has a 
mass of approximately 200kg. The following question pertains to the 
inertial forces experienced by this spinning rigid body, with I^y = Ixz = 

lyZ Xg Zg 0. 

Because of high rotational speed, the turbine is exceedingly sensitive to 
lateral balancing. Estimate the net force that the radial bearings would 
have to support for a 1mm mass imbalance, say, Pg = 1mm with respect 
to the spin axis. 

2. Kinematics. Although dead-reckoning is not recommended as a method 
for long-term navigation, it can be quite useful for pre-planning trajec- 
tories. For example, consider a surface vessel heading prohle as follows: 



t = 0 


0(f) = 0 rad 


0 < t <= 60s 


0(f) = t/00 X 7t/4 rad 


60s < t < 240s 


4>{t) = 7t/4 rad 


240s < t < 300s 


X 

o 

1 

o 

o 

CO 


t = 300s 


0(f) = 0 rad 





What is the total distance traveled in the Cartesian x direction during 
these five minutes? You may assume that the forward speed of the vessel 
is u = Im/s, and that there is no signihcant sideslip. Hint: Break the 
trajectory up into three parts, which can be evaluated separately. 





164 



24 PROBLEMS 



3. Coefficients and system modeling. Consider a weather vane in a 
wind of velocity Uo- If 0 is the angle of the vane with respect to the wind 
direction, 

(a) Write the single-degree of freedom [N) linearized equations of mo- 
tion about the fixed axis 0 . 

(b) Write Nq, Nq, and Ng in terms of N^, N^, W, etc.. 

(c) If we consider the differential equation 

Ay{t) + By{t) + Cy{t) = 0, 

the condition for stability is that A, B, and C must have the same 
sign. Express this requirement in terms of the derivatives in the 
previous question. Give physical interpretations for what would 
make such a device stable or unstable. 

(d) Create a numerical model of this system, using the MATLAB ODE 
solver ode45. The system equation can be written as two first-order 
equations: 



d^i 9 
dt [ 9 



-B/A -C/A 
1 0 




Simulate the system response to nonzero initial conditions (e.g., 
9(0) = 1,9(0) = 0). Discuss, using several examples, response 
sensitivity to B and C, which are related to the aerodynamic coef- 
ficients. Eor example, look at the range {A, B, C} = {1, ±3, ±3}. 




u 



0 
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4. Hydrodynamics. The figure below shows some characteristic fiuid 
force curves versus a motion parameter. Give the linear hydrodynamic 
coefficient at two different operating conditions, origin O and A: is it 
zero, small, finite positive, finite negative? 





5. Second-order system. Consider the vibration of a mass-spring-dashpot 
system: 



mx + bx + kx = 0, 

wherein all the coefficients are positive and real (i.e., physical). 

(a) Set x(t) = and find two solutions for s in terms of m, b, and k, 
using the quadratic formula. The solution pairs depend on whether 
b exceeds a critical value: what is the critical value and how do the 
solutions change? 

(b) From the initial condition a;(0) = Xo, and h(0) = 0, find the coeffi- 
cients for the general solution 

x{t) = + C26^'^. 

in terms of Xo, <si, and S2- What do you think will happen when 

Si = 52? 

(c) Write this x{t) for sub-critical b (that is, complex s) in terms of the 
standard parameters 
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CJn 

C 



(^d 

Give your answer in terms of a real exponential multiplied by a 
single sine or cosine, and make a sketch. 

(d) Sketch the response x(t) for supercritical b, i.e., both s are real, for 
the same initial conditions. What happens when the roots Si are 
far apart vs. when they are close together? 

(e) Consider now the case where an input acts to drive the mass from 
zero initial conditions: 



m 

b 



(undamped natural frequency) 
(damping ratio) 



2\Jkm 

oJn\J'^ — (damped natural frequency). 



mx + bx + kx = u. 

This equation dehnes a system, with input u, and output x. Using 
the Laplace transform, write the transfer function for this system, 
i.e., the impulse response, both in frequency (Laplace) space, and 
in the time-domain. Make a sketch of the time-domain result. Be 
sure to cover both the sub-critical and super-critical damping cases. 

(f) What is the step response of this system, in both the subcritical 
and super-critical damping cases? Include sketches. 

6. Submarine roll dynamics and control. A submarine has weight 
1200t (tons) and the center of gravity is 0.5m above the center of buoy- 
ancy (What can you conclude?). The rolling motion can be assumed 
to be decoupled from the other motions. This submarine has anti- 
rolling fins to ensure stability. The control hydrodynamic derivative is 
Ks = —2.8tm per degree of £n rotation 6, at a forward speed of 5m/ s. 

(a) Write the equation of motion for roll (K). 

(b) If the automatic control law S = ki'tp is used, where ip is the roll 
angle, what range of ki ensures stability? 
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(c) If the speed suddenly drops to ‘l.hmj how does the range of sta- 
bilizing k\ change? 

7. Steering. Draw curves of rudder angle S vs. yaw rate r for a dynamically 
stable surface vessel, and then for an unstable vessel. Indicate areas 
where the vessel turns against the rudder action. 

8. Submarine pitch/heave dynamics. Consider a submarine moving 
forward at speed U, and restricted to small motions in the vertical {x — z) 
plane. Assume that: 

• The submarine is symmetric in the x — y plane and yc = 0. The 
submarine is not symmetric in the x — z plane due to the sail, etc., 
but assume it is nearly symmetric, so that you can omit certain 
hydrodynamic terms. 

• At rest, the submarine is neutrally buoyant and stable with the x 
and y axes horizontal. What does this tell us about the magnitude 
of the buoyant force and where it acts with respect to the center of 
weight? 

(a) Derive the equations of motion for the submarine moving at speed 

U. 

(b) Derive the hydrostatic, restoring pitch moment for small pitch 6. 

(c) Linearize the inertial terms in the equations of motion. 

(d) Expand the fluid forces and moments in terms of the motions. Omit 
nonlinear and memory effects, and be sure to include the hydrostatic 
moment. Explain your choices. 

(e) Write out the complete linearized equations of motion. Does surge 
decouple from pitch and heave? Write out the coupled equations in 
matrix form. 

(f) Can the submarine be stable in pitch without feedback control? 
Can the submarine be stable in depth without feedback control? 

9. Submarine stability via slender body theory. 

(a) An underwater vehicle hull has the shape of a circular cylinder, 
except at the very ends, where it rapidly tapers down to a point. If 
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the radius is 0.5m, length is 17m, and speed is 2m/s, use slender- 
body theory to express the following coefficients: Mg, and M^, 

all referenced to the body midpoint. 

(b) It turns out that in order to match the lift and moment of the hull 
in experiments, a flat tail of radius 0.4m has to be included in the 
slender-body estimate. Determine the corrected value for M^,, and 
in addition compute Where is the aerodynamic center with 
respect to the center of the hull, and is this hull stable in pitch? 
Recall that the aerodynamic center is the equivalent location of the 
observed lift force that creates the observed moment. 

(c) A £n with pre-determined lift slope of 3.6/rad is to be applied at 
the tail to make the device just stable in pitch. How big is the 
required fin area that brings the net to zero, as desired? 

10. Lifting surfaces. Using the lifting surface formulas, estimate the coef- 
ficient Ky of a submarine sail as shown below, for U = lOm/s. 




11. Slender body theory. Consider a long ellipsoidal body of length I and 
diameter d. 

(a) With l/d = 7.0, approximate the cross-body added mass, using 
slender-body theory, and compare it with the exact results from 
the table below (Blevins). 

(b) Perform the slender body calculation also for a sphere, and compare 
again to the exact result. 

(c) What is the added mass in the in-line direction? 
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Sphere added mass 




2p7ia^/3 


Ellipsoid cross-body added mass 




Aapnab"^ / 3 


where a/h = 0.1 
0.2 
0.6 
1.0 
2.0 

5.0 

7.0 
10.0 

oo 




a = 0.075 
0.143 
0.355 
0.500 
0.704 
0.894 
0.933 
0.960 
1.000 


Ellipsoid longitudinal added mass 




Aapnab‘^/3 


where a/h = 0.1 
0.2 
0.6 
1.0 
2.0 

5.0 

7.0 
10.0 

oo 




a = 6.148 
3.008 
0.908 
0.500 
0.210 
0.059 
0.036 
0.021 
0.000 




12. Submarine pitch stability with various methods. We will calculate 
the pitch derivative due to an angle of attack, M^,, for a small submarine, 
using several different methods. 

The submarine is a rotationally symmetric ellipsoid with L/D = 7, L = 
35m. We £x the body origin to the midpoint as shown below. The body 
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is appended with two fins as shown in the figure: each has an area of 
geometric aspect ratio of 3, and has its longitudinal center of action 
Im forward of the body stern. 



(a) Apply wing theory to characterize due to the fins alone. 

(b) Apply slender body theory to estimate for the body alone. 

(c) Linearize the Munk moment, and give a corrected slender-body 
value for of the body. 

(d) Apply Jorgensen’s approximate formulas to estimate M^, for the 
body; make linearizations where necessary. 

(e) Use the experimental data from Hoerner (p. 13.2) below to estimate 

for the body. 

(f) Use the Hoerner result and your result from wing theory to write a 
net Mu, for the body plus fins. 

(g) What is the location of the net aerodynamic center? If it has an 
unstable location forward of the midpoint, what increase in fin area 
would bring it back to the midpoint? 

(h) What is the diameter of a fiat stern that will allow the slender body 
theory to give the same moment as the Hoerner (experimental) 
data, for the body alone? 

(i) What is the diameter of a fiat stern that allows the approximate 
formulas of Jorgensen to give the same moment as the Hoerner data, 
for the body alone? 



Hoerner p. 13.2; 

Symmetric body of revolution with L/ D = 6.7. 

Force and moment referenced to body midpoint: 

Z = —Q.bpU'^ D‘^Cydb^8-^{W /U) . Cy = 1.20. 

M = 0.5pU2DMC'„fetan(M;^). = 0.53. DESTABILIZING 
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13. Cable dynamics. A towing cable is steel-jacketed, with diameter 3cm, 
an effective extensibility of E = 100 x lO^Pa, and a density of pc = 
5000kg/m^. The vehicle it tows is streamlined, and has negligible mass 
(M). If we are towing at low speed and L = 4000m depth, 

(a) Calculate oji, the first undamped natural frequency in heave, with 
mL » M . 

(b) If the vessel heaves with a P = 2m amplitude and a 4.5-second 
period, what is the dynamic tension amplitude at the top (s = L) of 
the cable? The formula for this dynamic tension (with mL » M) 
is T = EAPksin{ks)/cos{kL), where A is the cross-sectional area 
of the cable, k = oo^m/EA, and m is the cable mass per unit 
length. 

(c) Taking into account the static tension induced by in- water weight, 
does the cable unload at the surface, for the heaving conditions 
above? 

14. Cable mechanics (hard!). You are asked to assess the operational 
envelope of a cable/vehicle system which has been installed on a vessel. 
The cable is steel-jacketed, with diameter 3cm, an effective extensibility 
of 100 X lO^Pa, and a density of 5000kg/m. The vehicle is streamlined, 
and has a mass (material plus added) of 100kg. 

(a) The towing angle at the vessel cannot exceed 25deg from vertical, 
for reasons of deck and crane geometry. If this angle is considered 
equal to the critical angle 0c, what is the maximum speed for the 
system? Assume the normal drag coefficient of the cable is Cn = 

1 . 2 . 
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(b) Considering the undamped, uncoupled axial dynamics, express the 
natural frequency as a function of the vehicle depth, and sketch the 
curve. If the fastest waves that the surface vessel responds to have 
period 4s, what is the “threshold” operating depth? 

(c) The undamped, uncoupled lateral dynamics follow the equation 



(m + rria) 



(fq 

dP 



^l)s 



+ Wn sin 



0 



dq 

(9s’ 



where rria is the added mass of the cable per unit length, and Wn is 
the in-water weight of the cable. Separation of variables g(s, t) = 
g(s) cos ut gives 



(9s2 



+ Wn sin 




+ Wnjq 



0, 



where 7 = (m -|- ma)oj‘^ /wn- Simplify this expression by writing 
T = WnS (the vehicle weight effect is small), sinc^ — 1 (the cable is 
nearly vertical), and then use the substitution 2 : = 2 ^/q^ to arrive 
at the Bessel equation 



2 ^ 

dz"^ 




+ z^q 



0 . 



The derivatives can be transformed from s to coordinates via the 
chain rule: 



dq 


dq 


dz 








ds 


dz 


ds 










d^q 


(dz^ 


1 dz 


d^z 


dz 


ds^ 


dz"^ 


[ ds j 


dzds 


ds 



(d) The equation above has a solution of the form q = cJq{z), where 
c is a constant to be found, and Jq{z) is the zero’th order Bessel 
function of the hrst kind. If q{s = L) = Q, that is, we impose a 
harmonic input at the top, hnd c in terms of Q and z{s = L). What 
condition makes q blow up, indicating a resonance frequency 
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(e) Recalling that strumming occurs when ujn — 0.2U/d, construct a 
graph of towing velocities vs. deployment depths for which the 
strumming would be centered. Include at least the three lowest 
modes in your sketch, and discuss the trends for the higher modes. 
Jq{x) = 0 for for x ~ 7r(n — 0.25), n = 0, 1, 2, • • •, especially when x 
is large. You should hnd that strumming is hard to avoid in deep 
towing applications! 



15. Vessel heading control. We consider the yaw/sway dynamics of a 
high-speed container ship. All the calculations below are to be made in 
nondimensional coordinates, so you do not need to perform any trans- 
formations in this question. 



(a) The nondimensional system with states x' = evolves accord- 

ing to dx' / dt' = Ax' + BS, where 



A 



-0.90 -0.42 
-4.8 -2.3 ’ 



-0.13 

1.4 



If the output is yaw rate r', i.e., C = [0 1], write the transfer 

function r'(s)/5(s) = C(sl — A)~^B. Use matrix inversion for the 
2x2 matrix (si — A), perform the matrix multiplications, and then 
use the quadratic formula if necessary to express your answer as 



r'(-s) _ S + Zi 

d(s) (s + pi)(s + P 2 )’ 

where 77 is a constant gain, —zi is the zero (there is one real zero), 
and [— pi, — P 2 ] are the two poles of the second-order system. Are 
both poles stable? 

(b) For the purposes of autopilot design, the transfer function (/>(s)/d(s) 
is crucial. Fortunately, all you have to do is divide r'(s)/<5(s) by an- 
other factor of s to account for the fact that d<p/dt' = r' . Sketch all 
the poles and zeros of the plant transfer function P(s) = 0(s)/<5(s) 
on a root locus plot. 




174 



24 PROBLEMS 



(c) A PID controller C'(s) will work to stabilize the heading angle, 
and brings to P(s)C(s) one additional pole at the origin, and two 
arbitrary zeros. Demonstrate a rough PID design by devising two 
controller zero locations that will attract the closed-loop poles into 
the left-half plane. Draw the path that you expect each closed-loop 
pole to follow, as the control gain increases from zero to very large 
values. 

You do not need to calculate the PID gains that go with your pro- 
posed zero locations. 

16. Full-state feedback of a submarine. The coefficients governing the 
pitch/heave dynamics of a Deep Submergence Rescue Vehicle (DSRV) 
are given below. They are nondimensionalized as in the lecture notes, 
using the factors p/2, U, and L. 



= 0.000118 
= 0.00193 

Xq = 0 

m' = 0.0364 
U = 2.0m/s 
L = 15.0m 



M'^ = -0.0113 
M'- = -0.00157 
= 0.0112 
M4 = -0.000146 
M's = -0.0128 
M'q = -0.156/7/2 



Z' = -0.0175 
Z'. = -0.000130 
z; = -0.0439 
Z/ = -0.0315 
Z's = -0.0277 



(a) Write the linearized (nondimensional) dynamics in the matrix form 
X = Ax + Bu, where the input u = 6, and the state vector is 
X = [w' ,q' ,6, Z']. Z' here is the elevation of the vehicle in inertial 
coordinates; your approximation for Z' should include both w' and 

e. 

(b) Is the DSRV vehicle open-loop stable, that is, without control action 
(5? You can assess this either by hnding the eigenvalues of your A- 
matrix above, or by computing the stability parameter C. 

(c) In preparation for controller design, create a simulation using Mat- 
lab. Make a graph of the open-loop step response, showing all the 
dimensional state variables versus dimensional time. 

(d) Assuming that all the states can be measured accurately, a full- 
state controller S = —Kx can be used, where 77 is a 1 x 4 gain 
matrix. Under this control, the system dynamics are governed by 
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X = [A — BK)x, so that li A — BK has all negative eigenvalues, 
the system is stable. 

Use the Matlab function place to put the four closed-loop poles at 
5 f(— 0.95 ± 0.31f) and 5 f(— 0.59 ± O.Slf), g = 1- Demonstrate that 
your closed-loop design is stable against nonzero initial conditions. 
What are the effects of increasing or decreasing g? 

Note that the pole locations suggested above are with respect to 
the dimensional system, i.e., using dimensional time. 

17. Nyquist stability. The inverted pendulum shown below is often used 
as a simple model for rocket flight, and can also illustrate the dynamic 
behavior of an unstable ocean vessel which is propelled from the stern, 
e.g., a barge being pushed by a tugboat. 

For this problem, we assume that all of the mass (m) is concentrated 
at the distal end of an arm of length /; a high-performance positioning 
system sets the horizontal position of the cart Ug- 




A Lagrangian derivation of the dynamics gives: 

0 = y sin(0) - y cos(0)iio. 

Note that the acceleration of the cart is now considered to be the input 
to the plant, e.g., u = iio- 

(a) If the observation is taken to be the angle of the bar 0, i.e., y = 
0, write the state-space representation of the dynamics linearized 
about the point 0 = 0. 
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(b) Show that this plant is fully state observable and controllable. 

(c) Now make the assumption that / = 1, and = 1 (it can be achieved 
through a proper nondimensionalization) What is the transfer func- 
tion of the plant, P(s)7 Use this plant model for the rest of this 
question. 

(d) Where are the plant poles and zeros, if any? 

(e) A stabilizing controller is not hard to hnd, and a suggested one is: 

cw = zlA±i), 

^ ^ s + 2 

Create a rough Nyquist plot of the loci of P(s)C(s), for the frequen- 
cies uj = [0, 2.2, cx)]ra(i/s. These three frequencies are all that you 
need to sketch the overall shape; there are no hidden loops or other 
features. In the case of u = 2.2, for which you will need to make 
explicit calculations, you may hnd the following identity useful; 

a + jb {a + jb){c — jd) {ac + bd) + j {be — ad) 
c + jd {c + jd){c — jd) 

Be sure to include the complex-conjugate points for —u on your 
plot. 

(f) Invoke the Nyquist stability criterion to conhrm that the closed- 
loop system is stable. 

(g) About how many degrees of phase margin does this design provide? 
What reduction in low-frequency gain can be tolerated? Is it a 
reasonable design? 

18. Root locus and loopshaping. The parameters governing the surface 
maneuvering of a high-speed container ship are given below for reference: 



m' 


0.00792 


r 

zz 


0.000456 


x'a 


-0.05 


L 


175.0m 


u 


8m/ s 


W 

V 


-0.00705 


y; 


0.0000 


Y' 


-0.0116 


Y' 

-l y. 


0.00242 


n 


-0.00258 


N' 


0.0000 


N/ 


-0.000419 


K 


-0.00385 


N' 

r 


-0.00222 


N's 


0.00126 
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Note that the center of vessel mass is located aft of the origin; for this 
model, the origin coincides with the center of added mass, so that 17 = 
Ni, = 0. The nondimensional system with states x' = [n7r'] evolves 
according to dx'/dt' = Ax' + B5, where 



A 



-0.90 -0.42 
-4.8 -2.3 ’ 



-0.13 

1.4 



The relevant output is yaw rate: C = [0 1] and D = 0. For the purposes 
of autopilot design, however, the transfer function 0(s)/<5(s) is needed. 



The following steps create two heading autopilots, using the root-locus 
and loopshaping techniques. In addition to the Matlab commands listed 
below, you will hnd very useful the convolution function convO which 
can be used to combine systems, e.g., for numerators, 
numPC = convCnumP ,numC) ;. Also, be sure that you equalize axis scal- 
ing for your plots in the complex-plane, by using axis ( ’ equal 0 ; . 



(a) Use the Matlab command tf (), or ss(), to create a system model 
of the open-loop transfer function P(s)C(s), using the plant above 
and a PID-type controller: 



C*(s) — kp(^ + Td-s + • 

The actual numerical values for kp, Td, and are to be found in the 
next step. 

(b) Using Td = 2 and Tj = 6 as suggested values, use the Matlab com- 
mand rlocusO and then rlocfindO to select a controller gain 
kp, that puts the three slow poles in the following sector: 1) min- 
imum undamped frequency (nondimensional) of 0.3, 2) maximum 
frequency of 0.5, and 3) minimum damping ratio 0.7. Give a root 
locus plot, with your pole locations clearly marked on top of the 
trajectories taken as kp varies. You don’t need to show the fourth, 
fast pole, which will be quite far to the left. 

(c) Apply the kp you selected to P{s)C{s), and then use the Matlab 
functions feedback () to create the resulting feedback system, and 
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stepO to plot the closed-loop system response to a step input in 
desired heading. 

(d) Use the Matlab command nyquistO to make a Nyquist plot of 
P(s)C(s) for your design. Make a visual estimate of the gain and 
phase margins. 

(e) An alternate approach for controller design of this stable plant 
is loopshaping: For the open-loop function L(s) = uic/s, where 
ojc = 2.0, invert the plant to come up with a compensator: C(s) = 
L{s)/P{s). This design has inhnite gain margin and 90 degrees 
phase margin. 

(f) As above, create the feedback system, and plot the closed-loop step 
response. 

(g) The loopshaping control is not quite a PID-controller; how does it 
differ, and what would L{s) have to contain to make it a PID? 

(h) The controllers you just designed are in nondimensional time coor- 
dinates; give the P,I, and D gains for use on a real time scale, for 
the root-locus design. 

19. LQR. Consider the state-space system and LQR design: 



B = [1 0]^ 

C = [0 1] 

D = 0 
Q = C^C 
R = p. 

The plant is an undamped oscillator with undamped poles at ±j. Note 
that the plant output is position for this problem. 



What is the control gain K in terms of p? Hints: There are two 
solutions for p^', choose the positive one. Also, the expression for 
P 22 is messy; luckily, you won’t need to use it. 
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(b) Determine the limiting approximations for K with p very small 
and very large - these are the cheap control and expensive control 
problems. 

(c) Derive the limiting closed- loop pole locations for p — 0, giving the 

frequency and damping ratio of the Butterworth pattern in terms of 
p. You can get the characteristic equation for the poles as det{sl — 
{A — BK)) = 0, and then make it £t the form -|- 2C,ojn -|- = 0. 

20. LQR. Consider the first-order, unstable system governed hy x = x + u, 
for which the output is y = x. The state-space matrices are A = B = 
C = 1. 

(a) Solve the LQR Riccati equation for P, for the case Q = C^C and 
R = p. Your answer should give a positive P as a function of p 
alone. 

(b) Consider the cheap-control problem, with p — 0. Write the leading- 
term approximations for P and then K, and compute the eigenvalue 
(there is only one) of the closed-loop system. Note that eig{M) = 
M, if M is a scalar instead of a matrix. 

(c) Now look at the expensive-control problem, with p ^ oo. Again, 
write the approximations for P and K, and find the closed-loop 
system eigenvalue. 

21. LQG/LTR. For the inverted pendulum plant (see previous question) 
with state-space matrices 



A 



0 1 
1 0 ’ 



C = [0 1] , 



the KF Riccati equation yields H = [100 14]^, for the choices Vi = 
[100 0 ; 0 0], and 1/2 = 0.01. 



(a) Compute the open-loop transfer function for the Kalman filter loop: 
L{s) =C{sI - A)-^H. 
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(b) Even though this L(s) is unstable, its magnitude plot conhrms that 
it is a reasonable loopshaping design. What is the low-frequency 
gain for this L(s), and what kind of tracking performance should 
we expect from the associated LTR design? 

(c) Now consider the closed-loop transfer function, which you can write 
as S(s) = L(s)/(1 -|- L(s)) (the sensitivity). What is the character- 
istic equation for the closed-loop system, and about what damping 
ratio has the Kalman hlter provided? 

22. LQG/LTR design. The parameters for the linearized sway/yaw mo- 
tions of a swimmer delivery vehicle are given below. 

% Parameters, all nondimensional except [U,L] 



u = 


4.0 ; 7o m/s 


L = 


5.3 ; 7„ m 


Izz = 


0.006326 ; 


m = 


0.1415 ; 


xg = 


0. ; 


Yv = 


-0.1 ; 


Yr = 


0.03 ; 


Nv = 


-0.0074 ; 


Nr = 


-0.016 ; 


Ydelta = 


0.027 ; 


Yvdot = 


-0.055 ; 


Yrdot = 


0. ; 


Nvdot = 


0. ; 


Nrdot = 


-0.0034 ; 


Ndelta = 


-0.013 : 



You are asked to develop an LQG/LTR controller for this plant, and it is 
suggested that you compose a single Matlab script to perform the steps 
in sequence. Please make sure you answer all the questions, and include 
a listing of your code. This entire design is made in nondimensional 
coordinates. 



(a) Plant Modeling and Characteristics 
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i. Construct a state-space plant model, to take rudder angle 6 as 
an input and give heading angle 0 as an output. Please provide 
the numerical values for the A, B, C matrices. There should be 
three states in your model, with one input channel and one 
output channel. 

ii. Compute and list the controllability and observability matrices; 
is the plant state-controllable and state-observable? 

hi. Where are the poles of your plant model? Is this model stable? 
iv. Show a plot of your plant’s step response. 

(b) LQR and KF Designs 

i. Using the Matlab command lqr(), you can compute the LQR 
feedback gain K, for given A, B, Q, and R matrices. With the 
choices Q = C^C, and R = p, list K and plot the closed-loop 
step responses for the choices p = [0.1,0.001,0.00001]. How do 
the gains and step responses change as you make p smaller and 
smaller? 

Note that the fundamental closed-loop LQR system is 



X BK^x BKxd^giy-^d 

y = X, 

i.e., the input to the closed-loop system is x desired and the out- 
put is X. Your plot should show specihcally the output 0, for an 
input of Xdesired ['f^desired ^desired ^■/(pdesired 1] • This 
compression can be achieved in one step by premultiplying the 
system by and post-multiplying it by C: 



X = (A - BK)X + BKC^Hdesired 

y = Cx, 

ii. The Matlab command lqe() can be used to generate the Kal- 
man Liter gain H, given design matrices A, C, Vi, and V 2 . For 
the choices Vi = Jsxs and V 2 = 0.01, compute H, and make 
a plot of the closed-loop step response. Be sure to give the 
numerical values of H. 
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Note that the lqe() command asks for a disturbance gain ma- 
trix G; you should set this to /sxs- The closed-loop KF system 
is as follows: 



k = {A-HC)x + Hy 

y = Cx, 

i.e., the input is the measurement y and the output is an esti- 
mated version of it, y. 

(c) Loop Transfer Recovery 

The LQG compensator is a combination of the KF and LQR designs 
above. With normal negative feedback, the compensator C'(s) has 
the following state space representation: 



T = {A- BK - HC)z + He 
u = Kz, 

so that the input to the compensator is the tracking error e = 
r — y, and its output u is the control action to be applied as 
input to the plant. The total open-loop transfer function is the 
P{s)C{s)] in Matlab, you may simply multiply the systems, e.g., 
sysPC = sysP * sysC ;. 

i. Make a log(magnitude) plot of the KF open-loop transfer func- 
tion L{s) = C{sl — A)~^H, versus log(frequency). You may 
hnd the Matlab command freqrespO helpful. |L(s)| should 
be large at low frequencies, and small at high frequencies, con- 
sistent with the rules of loopshaping. 

ii. As p — 0, the product P{s)C{s) L(s). Demonstrate this by 
computing P{s)C{s) for the three different values of p above, 
and overlaying the respective |P(s)G(s)| over the plot of part 
3a). 

hi. Make a closed-loop step response plot for the smallest value of 
p. How does it compare with the KF step response of part 2b)? 

In real LTR applications, the particular values of Vi and V 2 can be picked 
to control the low-frequency gain, and crossover frequency of the open- 
loop KF system L(s) = C{sl — A)~^H. 
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23. Inertial navigation. A triaxial accelerometer package is aligned with 
the [x,y,z] (fwd, port, up) axes of an underwater vehicle. At a particular 
instant in time, the three raw accelerations from the strain guages are 
[x = — 0.96,i/ = 0.80, i = 9.73] meters/second. Note that under the 
influence of gravity alone, this sensor will always report that the vehicle 
is accelerating upwards. 

(a) Ignoring the coupling of Euler small angles, what estimates can 
you give for the roll and pitch angles? I suggest that you draw the 
measurement vector in the coordinates of the vehicle frame, and 
then consider the orientation of the body which would lead to it. 

(b) Given that the acceleration due to gravity is 9.81m/s, how do we 
know that these are reasonable (but not foolproof!) angle calcula- 
tions to make? 
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24. LQG/LTR design: term project. You are asked to carry out mod- 
eling and control system design for the surge dynamics of an oceanic 
cable-laying vessel. This vessel needs to have very good speed control 
so that the cable, which is being paid out constantly and sinks slowly 
under its own weight, lies properly on the seafloor. The bathymetry 
for the path is well-known but variable, so the vessel will need to slow 
down and speed up frequently. Unexpected deviations in speed from the 
calculated trajectory are likely to create loops or kinks in the cable, or 
induce large tensions. 



vessel 


length 


L 


145m 




draft 


h 


9m 




beam 


b 


22m 




displacement 


V 


17750m3 




surge added mass 


ma 


0.07pV 




thrust reduction factor 


t 


0.19 




wetted surface area 




5585m^ 




towed resistance coefficient 


Cr 


0.0025 




wake fraction at propeller 


w 


0.22 


propeller 


diameter 


D 


7.21m 




pitch/diameter ratio 


P/D 


1.1 




rotative efficiency 


Vr 


1.025 




zero-speed thrust coefficient 


kt{J = 0) 


0.6 




slope of kt with J 




-0.522 




zero-speed torque coefficient 


kq{J = 0) 


0.1 




slope of kq with J 




-0.0833 




rotational moment of inertia 


Ip 


2.03 X lO^kgm'^ 


engine 


parameter 


a 


0.876 






b 


0.208 






c 


2.173 






d 


-.0939 




maximum fuel rate 


fm 


1.5kg/s 




maximum torque 


Qm 


70000Nm 




maximum speed 


'^m 


60Hz 


gearbox 


reduction ratio 


A 


32 




efficiency 


% 


0.97 
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It is desired to create a linear controller for the nominal surge motion, 
and then demonstrate that it will work with a simulation of the real 
nonlinear system, within a neighborhood of the nominal speed. 

The vessel is powered by a single gas turbine engine, driving one pro- 
peller. The parameters for the vessel, propeller, and engine system are 
as given in the table. 

The gas turbine torque-speed characteristic hts the relation: 










where Qe represents the engine torque in iVm, / the fuel rate in kg/s, 
and He the rotation rate of the engine in Hz. 

(a) Make a map of the gas turbine characteristics. For instance, make 
a contour plot which gives Qe/Qm (the fraction of maximum torque 
developed) as a function of f / fm (the fraction of the maximum fuel 
rate) and Ue/ Um (the fraction of the maximum rotation rate of the 
engine) . 

Be sure to include clipping on your contours, where the calculated 
torque exceeds Qm- Similarly, the engine cannot develop negative 
torque. 

(b) Make a table of some steady operating conditions for fo/ fm in the 
range of 0.05-0.95, in which you show: advance ratio seen by the 
propeller, engine torque, propeller rotational speed, engine rota- 
tional speed, and vessel linear speed. 

Plot the pairs [fo/ fm, neo/nm] on the characteristic plot from Part 
1, and assess whether the vessel, engine, and propeller are well 
matched to operate over a range of fo- 

(c) Construct a linear approximation for the plant dynamics around 
the operating condition fo/ fm = 0.8. If the plant input is <5/, and 
the output is bu, list the A,B,C matrices for your plant, and the 
eigenvalues. 5 here indicates the perturbed value from the steady 
state, e.g., f = fo + bf. 
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(d) Create an LTR-type controller for the plant. First, choose the hlter 
gains Vi and V 2 to make the crossover frequency (where the magni- 
tude of L{s) = C{sl —A)~^H passes through unity) about 0.6rad/s. 
This choice corresponds with the closed-loop system being able to 
overcome waves with period lOseconds or longer. Note that for this 
scalar design, you will be able to move the curve |T(s)| up or down 
on the Bode plot, but you cannot easily change its shape, or move 
it side-to-side. 

Secondly, recover this loop shape with an LQG-type controller, us- 
ing a small control penalty. 

You should prepare for this part: a plot of log{\L{s)\) versus log{u), 
a listing of your Vi and V 2 choices, and the A, B, C matrices of your 
LQG compensator. 

(e) In the event of a current or wind disturbance, the actual vessel 
speed will vary from the nominal value. This is not acceptable in 
the long term, given that the vessel is laying out a cable. 

An integral action can be added quite easily in the LTR design 
technique. First, add an integrator to the input of the plant, so 
that an “augmented plant” is created. This augmented plant has 
one more state than the original plant: 



Bang = [O 1]^ 

Caug = [C 0 ] 

The idea is to then carry out the KF design and LTR as before, with 
the augmented plant model. When this is done, move the integrator 
from the augmented plant into the compensator. The addition of an 
integrator channel to the plant or compensator can be accomplished 
with the command sysPaug = sysP * tf([0 1] , [1 0]) ;, forex- 
ample. 

Present a hgure of |T(s)| and the other quantities as you gave them 
in Part 4. This is to be kept as a separate controller design. 

(f) For the purpose of demonstrating your controllers, construct a non- 
linear simulation of the plant system, with input /. You don’t 
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need to show anything for this task; it is mostly to verify that 
the steady conditions you calculated above, with / = 0.8/m, are 
actually steady conditions from the point of view of the whole sim- 
ulation. You may wish to use global variables in Matlab so that 
you don’t have to type everything in more than once. 

(g) Augment your simulation above with the compensator dynamics. If 
the compensator has states and your simulation has states [m, np], 
make a new state vector [m, np, z\. You will propagate the first two 
states as usual with the nonlinear equations and the fuel rate / as 
input. The compensator states 2 : are propagated by a set of system 
matrices [Ac, 5c, C'c], which you just designed. Note that the input 
to the compensator is an error signal: e = r — u, where u is the 
speed of the vessel and r is a reference speed. The output of the 
compensator is then fuel rate /. 

We are interested mainly in disturbance rejection, for which you 
will just set r = Mq, the steady condition speed. Recall that you are 
implementing a controller for operation about a nominal condition; 
whenever your compensator generates / from the error signal, be 
sure to then add on 0.8/m, the part needed to keep us near the 
steady condition. 

Demonstrate four properties, showing time-domain simulations for 
each: 

• The hrst controller, without integral action, does not com- 
pletely reject steady disturbances. I suggest you implement 
a disturbance as an acceleration: 

uprime = uprime + k. In the absence of a controller, this 
line would drive the system with an additional acceleration of 
km/s^, thus modelling the effects of steady wind or current. 
This line is to be added after computing the usual parts of 
uprime. 

• The hrst controller does not reject high frequency disturbances. 
A suggestion is to let 

uprime = uprime + k*sin(10*omega_c*t), where omega_c is 
the crossover frequency. This adds an oscillating acceleration 
at ten times the crossover frequency. 

• The hrst controller does, however, reject some lower frequency 
disturbances. For example, try 
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uprime = uprime + k*sin(omega_c/10*t)/100. Note the ex- 
tra factor of 100 is included so that the result can be compared 
directly with the high-frequency disturbance case above. The 
effect of an acceleration disturbance onto the observed output, 
velocity u, scales with the inverse of the frequency. 

• The second controller, with integral action, rejects steady dis- 
turbances. Use the same approach as above for steady distur- 
bances 

What is the propeller response during the high-frequency distur- 
bances? Is it reasonable or do we need to slow down the controller 
bandwidth? 




